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For centuries, military forces have used camouflage to obscure potential targets 
from the enemy. Because the eye is fairly adept at picking out edges, colors, and 
bright areas, camouflage is often used to degrade these qualities from human detection. 
The purpose of this thesis was to investigate the role of certain spatial, temporal, and 
chromatic features on the human visual system and how these features may aid the 
quest for better camouflage. Methods: Test patterns were spatio-temporal raised 

cosines of varying orientation (horizontal or vertical and oblique), spatial frequency (1, 
3, and 7 cpd), and modulated at 2.0 Hz.- Color contrast thresholds were determined 
from 16 different red-green color mixture ratios. This methodology eliminates the 
problems with luminance artifacts and the need to determine exact equiluminance. 
Results: The data formed an ellipse with the half-length measuring color discrimination 
and the half-width measuring brightness discrimination. A maximum likelihood 
method was used to fit the data. Three of the four subjects showed a 3 cpd chromatic 
oblique effect, while the 1 and 7 cpd achromatic and chromatic oblique effect was 
inconsistent across subjects. Conclusions: While real-world objects are more complex 
than laboratory stimuli, knowledge of spatial and chromatic qualities that inhibit 
detection will aid the quest for better camouflage. 
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EXECUTIVE SUMMARY 



For centuries, military forces have used camouflage to obscure potential targets from 
the enemy. Because the eye is fairly adept at picking out edges, colors, and bright areas, these 
are the qualities that camouflage seeks to degrade. The purpose of this thesis was to 
investigate the effects of certain spatial, temporal, and chromatic features on the human visual 
system and how these features may aid the quest for better camouflage. 

Camouflage targets tend to match their backgrounds both in color and structure. 
Netting or paint can make detection of a potential target more difficult by reducing chromatic 
contrast. They can also reduce structural contrast by reducing sharp edge effects. This is 
evident in the "stealth" design in which the lack of defined edges helps reduce the radar 
signature. By removing defined edges, and thus high spatial frequency components, the visual 
signature is also reduced since "stealth" ships and aircraft generally have backgrounds 
consisting primarily of low frequencies (sea and sky). The perception aspects of this real world 
example can be simplified by examining less complex stimuli in the laboratory. For example, 
camouflage design would be affected if it were known that oblique chromatic lines were more 
difficult to detect than horizontal or vertical chromatic lines. 

It has been shown that horizontal and vertical achromatic lines are easier to detect than 
oblique achromatic lines. This phenomenon is known as the "oblique effect," which in this case 
is an achromatic oblique effect. Several studies have shown that the magnitude of the oblique 
effect is largest with high spatial and a low temporal frequency sinusoidal gratings. Previous 
researchers have used this knowledge to design experiments testing for a chromatic oblique 
effect, but have had problems with luminance artifacts due to the difficulty of obtaining exact 
equiluminance. By adapting the methodology of an earlier researcher the problem of 
determining exact equiluminance was avoided and an experiment to test for a chromatic 
oblique effect was designed. 

The experiment was conducted concurrently at the Naval Postgraduate School (NPS) 
and the University of Louisville, Kentucky (UL). Four subjects, 2 NPS and 2 UL, volunteered 
for this experiment. All subjects had normal (20/20), or corrected to normal, acuity and color 
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vision. Stimuli were presented by a VisionWorks computer graphics system (Vision Research 
Graphics, Inc.) on an IDEK MF-8521 high resolution color monitor (21" X 20" of viewable 
area). The monitor had a resolution of 800 by 600 pixels (x=75.02 and y=74.92 pixels/degree), 
98.9 Hz frame-rate, mean chromaticity of r = 0.334, g = 0.336, b = 0.300 (1931 CIE), and a 
maximum luminance of 100 cd/m2. The University of Louisville's apparatus and procedure 
were identical to the Naval Postgraduate School's, except that the stimuli were displayed on a 
17" Nanao Flexscan F2.21 color monitor. Subjects viewed the monitor from 1.5 meters and 
were positioned by an adjustable chinrest. 

Sinusoidal gratings were presented within a spatially windowed circular test field that 
subtended 7.59° of visual angle. The Gaussian window was truncated at ±1 standard 
deviations for both x and y directions. The test patterns were one-dimensional spatio-temporal 
sinusoids of varying orientation (principal and oblique), spatial frequency (1.0, 3.0, and 7.0 
cycles/degree), and color contrast. Test patterns for each subject consisted of two orientations, 
principal (0° and 90°) and oblique (45° and 135°). For each subject, maximum sensitivity for 
each orientation within the principal and oblique grouping was chosen. All sinusoids were 
raised cosines temporally modulated at 2.0 Hz. The sinusoid pattern was presented in a 1500 
msec interval with contrast ramped on and off according to a linear window. (Contrast 
peaked at 202 msec and fell at 1304 msec. Color contrast was computed by different ratios of 
percent red and green. Sixteen different sinusoidal red-green color mixtures were generated by 
changing the red phosphor only, green phosphor only, or by changing the red and green guns in 
fixed proportions. Color contrast was defined according to the Michel son formula. Thresholds 
were determined by a sequential two-altemative forced choice adaptive psychometric 
procedure, QUEST. Threshold was defined at 75 percent correct. A total of 480 trials, 30 
trials per condition, were randomly presented within each session. A session (~ 45 minutes) 
consisted of one sinusoidal condition with 16 different red-green color mixtures. A subject had 
to complete six sessions to contribute one set of 16 thresholds for each condition. 

Numerous surveys of differential thresholds have been carried out, but one of the more 
extensive ones was completed by D.L. MacAdam. The data from this survey was elliptical (the 
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closed curves connecting the thresholds were elliptical in shape). It was shown that the errors 
of thresholds about these closed curves were Normally distributed, therefore the curves should 
be ellipses, as they appeared to be. 

The elliptical properties of the experimental data were used to fit ellipses using the 
method of maximum likelihood a nested F-test (which is only approximate, because of non- 
linearity) was used to examine the significance of the different orientations. Results of this test 
showed that the chromatic oblique effect was inconsistent across the spatial frequencies of 1, 3 
and 7 cycles/degree (cpd). The lack of a 1 cpd chromatic oblique effect was due to the 
insensitivity of the chromatic channel at the lower spatial frequencies. Two of the four subjects 
showed a graphical 7 cpd chromatic oblique effect, but this was non-significant. Chromatic 
aberration may have been a factor. Three of the four subjects showed a 3 cpd chromatic 
oblique effect with two significant and the third marginally significant. It is predicted that the 
marginally significant subject would show significance with additional trials. 

The main value of this study is the tool it provides for further investigation of a 
chromatic oblique effect without the problems associated with luminance artifacts. 
Additionally, further investigation of a chromatic oblique effect will likely provide knowledge 
of spatial and chromatic qualities that inhibit detection. Knowledge of these qualities will aid 
the quest for better camouflage 
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I. INTRODUCTION 



The continual decline of the military budget has necessitated the increased protection of 
current and future war-fighting assets. This increase, coupled with public expectation of zero 
or close-to-zero casualties, has forced the services to reassess the way they conduct operations. 
Today's political climate necessitates this reassessment since the potential loss of public support 
for military actions generally increases with the number of American casualties. 

A cruise missile attack is an effective method used to minimize civilian and friendly 
casualties. For example, cruise missile attacks were successful against Iraq’s military targets 
during the Gulf War. Although the cruise missile is an effective weapon, it is an expensive 
resource for the United States military inventory. A less expensive alternative is an air strike, 
but the disadvantage of an air strike is the increased likelihood of aircrew casualties and missed 
targets. Ideally, the aviator would like to enter the threat zone undetected, thereby increasing 
the probability of locating the target without becoming engaged by the enemy. 

The ability to avoid detection is a distinct advantage during battle. For centuries, the 
element of surprise has resulted in a quick and decisive destruction of forces. For example, the 
U.S. Air Force F-117 “stealth” fighter was responsible for much of the precision bombing 
during the Gulf War. Because their aircraft was nearly invisible within the Iraqi air defenses, 
pilots had additional time to accurately drop bombs on target. This near-invisibility must 
extend beyond the visible spectrum since today’s battlefields are equipped with electro-optical 
sensors that often extend the range and increase the probability of detecting potential threats. 
Electro-optical sensors make detection possible through visual, infrared, thermal and other 
means. For example, shielding hot parts of a vehicle is part of the vehicle's thermal camouflage 
as it tries to disperse its thermal signature. 

But, even with these technological advances, a large threat to operations, and to 
reconnaissance operations in particular, is often the most common sensor~an enemy’s eyes. 
Prior to the start of the Gulf ground war, U. S. forces had reconnaissance teams operating in 
the interior of Kuwait. These teams used their speed and camouflage to prevent detection and 
capture. If their camouflage had been inadequate, capture and death would have been likely, 
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delaying the start of the ground offensive. While the saying “if they can't see you, they can't 
hurt you” is no longer true, the likelihood of being hurt is reduced when the enemy cannot see 
you. 

To make the enemy’s task of detection more difficult, camouflage is used. Since 
detection has traditionally been associated with visual detection, camouflage is generally 
thought of as making the visual detection of personnel or any potential target more difficult. It 
is in this sense that camouflage will be discussed. 

Camouflage has many different applications, ranging from the clothes and face paint of 
an infantryman to the netting used to cover tanks and vehicles to the paint used on larger 
platforms such as aircraft and ships. Since the eye is fairly adept at picking out edges and 
bright areas, camouflage is often used to break up edges and to cover or conceal bright areas. 
All objects possess certain unique qualities of shape and color. In order to deceive a sensor, 
the object must blend with the background. An invisible object would match the background, 
while an easily identifiable target would contain noticeable spatial, temporal, or chromatic 
features. By manipulating the spatial, temporal and chromatic features, an object can be made 
more difficult to detect. Knowledge of the range of these features will aid military designers in 
their quest for better camouflage. 

To understand how to make detection more difficult, one must also understand the 
sensors that will be used. Since, in this discussion, the primary sensors are the enemies’ eyes, 
either directly or through some sort of image intensifying mechanism, a working knowledge of 
how the eye works and what cues it uses to accomplish detection is essential. 
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II. BACKGROUND 



The ability to perceive an object is an unconscious, automatic process. The world is 
filled with a variety of sensory information that stimulates our senses. This information is 
obtained through visual, auditory, tactile and olfactory inputs. We use our senses to collect this 
information and translate it into meaningful units of sensory awareness. This information is 
then relayed to the brain, resulting in the formation of perceptions. The brain then categorizes 
this sensory data and compares it to past experiences. Thus, perceptions are a culmination of 
sensory inputs that are organized into a meaningful representation of the outside world. The 
study of perception involves a complete .understanding of the description of objects, 
appearances and events in the outside world (Sekuler and Blake, 1994). In brief, sensation and 
perception refer to a sequence of events: stimulation of an external object; machinery to 
capture this information; and translation of this information into electrical energy to form an 
experience. Perceptual experiences guide our actions in the world. This thesis investigates the 
different techniques of manipulating visual sensory input to change the appearance or 
perception of an object. (Sekuler and Blake, 1994) 

The basic function of the eye is to capture visual sensory input or light and focus it on 
the retina, a thin layer of receptor cells located in the back of the eye. The light must pass 
through these retinal cells before reaching the photoreceptors, which convert the light into an 
electrical signal. This process of converting an external stimulus into a neural signal is called 
transduction. These neural signals are then sent through the network of retinal cells, which, in 
turn, are sent to the brain (Figure 2.1). 

The human retina contains two major classes of photoreceptors, rods and cones. There 
are approximately 8 million cones and 120 million rods. The fovea, located at the center of the 
retina, has the greatest concentration of cones. The cones are sensitive to daylight and provide 
high-acuity color vision, while the rods are used for night viewing and are thought to be 
achromatic. Between five and ten percent of the total cone population are sensitive to short 
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Figure 2.1. A cross section of the human retina. From Sekuler and Blake [1994]. 
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wavelengths, with a peak sensitivity of 420nm. This cone class is referred to as S or blue 
cones. A few S cones are located at the fovea, but the largest concentration of S cones forms a 
ring around the fovea. This large concentration of cones tapers off with increasing distance 
(eccentricity) from the fovea. The remaining cone population is sensitive to middle (peak 
sensitivity of 530nm) and long wavelengths (peak sensitivity of 565nm). The long-wavelength 
(L or red) cones outnumber the middle-wavelength (M or green) cones two to one. The R and 
G cone distribution is randomly mixed throughout the retina, with the greatest concentration 
located inside the fovea. (Tovee, 1996) 

While he did not directly identify cones, Thomas Young believed that the retina 
contained three receivers that were sensitive to a limited number of light vibrations. These 
receivers are now known as cones. Young also is responsible for one of the first explanations 
of how we perceive color. Earlier, Newton had demonstrated that white light could be split by 
a prism into a spectrum of colored lights. He found that recombining some of these colored 
lights resulted in the original white light. Newton mixed and subtracted colors, but he did 
not attempt to explain how we perceive color. Young, however, postulated that color 
perception is due to the vibrations of light interacting with the retina. The three receivers in the 
retina were broadly tuned with overlapping sensitivities. Helmholtz confirmed and elaborated 
Young’s color theory by showing that there are three types of receivers (cones or 
photoreceptors) in the human retina, and that each type contains a different pigment. The 
spectral sensitivity of the cone is determined by the absorption spectrum of its photopigment. 
(McDwain, 1996) 

Young hypothesized that there were three broadly tuned receivers in the retina 
because, he reasoned, a single broadly-tuned receiver could not provide enough information 
about the wavelength of light. If there were two receivers, then there would be one particular 
frequency that excited both receivers equally, and thus white light would be produced 
(intersection of the curves in Figure 2.2b). Since Young did not observe white light in the 
color spectrum produced by a prism, he concluded that the visual system must have three 
broadly tuned receivers. A single cone pigment cannot discriminate between changes in 
wavelength and changes in the intensity of light. A cone can only increase or decrease its 
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output, so its signal is ambiguous as to whether the change is due to a shift in the light’s 
wavelength or to a change in its intensity. This type of response is explained by the principle of 
univariance. For example, a retina that contains one cone class may give the same response to 
two different wavelengths, while a two-cone class retina may excite the receptors in different 
ratios (Figure 2.2). Therefore, the retina needs at least two cone types to distinguish between 
changes in wavelengths. Primate vision can discriminate millions of colors with three cone 
types, whereas non-primate mammals, which usually have two cone types, cannot. These 
mammals tend to rely much more heavily upon auditory and olfactory sensory inputs to 
survive. (Tovee, 1996) 

A B 





Figure 2.2. Wavelength discrimination by (a) a one-class retina and (b) a two-class retina. 
A retina that contains one cone pigment responds more or less with the same energy for 
some wavelength within its spectral sensitivity~the principle of univariance. A retina that 
contains two cone pigments will have different responses depending upon the location of 
the two wavelengths. From McUwain [ 1996]. 

The Young-Helmholtz trichromatic theory accounts for many, but not all, of the 
phenomena associated with color vision. The theory predicts that a signal comprised of a 
certain combination of long and medium wavelengths cannot be distinguished from a specific 
third wavelength (yellow), but it does not account for the fact that the signal appears yellow 
rather then red-green or green-red. Additionally, the phenomenon where an object’s color 
appears to vary depending on the colors viewed immediately before viewing the object 
(successive color contrast) or on the colors surrounding the object (simultaneous color 
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contrast), cannot be accounted for. A theory that does account for both these phenomena and 
those explained by the Young-Helmoltz theory is known as the opponent-color theory. 
(McDwain, 1996) 

The opponent-color theory was first introduced in 1878 by Hering and has been 
furthered through the work of Hurvitch and Jameson. Hering postulated that there were three 
visual processes, two chromatic and one achromatic. The processes consist of three 
antagonistic or opponent pairings. These pairings are red-green and yellow-blue for the 
chromatic processes and black-white for the achromatic process. Such opponent pairs are well 
explained by the center-surround or on-center and off-center type ganglion cells consisting of 
the aforementioned pairings. These opponent pairs account for the fact that the colors in these 
pairings cannot be seen at the same time; thus, there are no reddish-green hues. The inputs 
from the S, M and L cones in the first diagram of Figure 2.3 display in a simplified way how 
these inputs are combined to form a signal that we perceive as blue. The second diagram in the 
Figure is identical to the first, except that the weightings of these signals result in the perception 
of yellow instead of blue. Intermediate weightings of the signals displayed result in our 
perception of many more colors than the six shown in Figure 2.3. (Tovee, 1996) 




“SUE" "YELLOW" "GREEN" ' "RED" "WHITE" "BLACK" 

Figure 2.3. Simplified hypothetical display for a model based on color opponency. From 
Mcllwain [1996], 

The cones relay their information by synapsing onto ganglion cells, whose axons travel 
through the optic nerves to the visual cortex located in the back of the brain. The human eye 
contains approximately one million ganglion cells. The input of 128 million photoreceptor 
signals has been reduced to an output produced by these one million ganglion cells. In 1938, 
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Hartline discovered that the retinal ganglion cells of the frog were comprised of two 
concentrically shaped ring-like areas on the retina (Figure 2.4). He found that certain ganglion 
cells increased their electrical energy when a light was passed through their center or inner ring 
and decreased their electrical energy when a light was passed through the outer ring. Hartline 
called these ON-center ganglion cells. Cells that respond vigorously to a dark center and a 
light edge are called OFF-center ganglion cells. A network of such ON-center and OFF-center 
cells is responsible for providing edge detection as well as orientation information. Hartline’ s 
work laid the foundation for later work by Devalois (1958), who discovered that primate 
ganglion cells behaved in a similar manner. (Sekuler and Blake, 1996) 

The morphological and physiological properties of primate ganglion cells are 
divided into three categories: large size P a or A cells, small size P p or B cells, and P-, or W- 
like cells. Livingstone and Hubei (1988) classify two major types of ganglion cells, P a and 
Pp cells, and their projections to different cortical regions within the primate visual system. 
P a cells (ten percent of ganglion cells) have high conduction velocities, transient 
responses, no color sensitivities, high contrast sensitivity and very good temporal 
frequency modulation. P p cells (80 percent of ganglion cells) have lower conduction 
velocities, sustained responses, low contrast sensitivity, color opponency and moderate 
temporal resolution. Both P a and Pp neurons are segregated into two different pathways 
which project to different locations within the lateral geniculate nucleus (LGN), VI, and 
higher cortical regions within the primate cortex (Livingstone and Hubei, 1984). Refer to 
Figure 2.5 for a graphical representation of the visual pathway. (Merigan, 1989) 

The magnocellular pathway (M pathway) receives input from P a ganglion cells that 
project first to layers 1-2 of the LGN. Cells in the magnocellular geniculate layers project 
to layer 4C a of the primary visual cortex. From layer 4C a they project to layer 4B, which, 
in turn, projects to visual area 2 and to the medial superior temporal area (Tovee, 1996). 
This pathway is thought to involve spatial awareness, that is, ‘where’ an object is located 
in space. Alternatively, the parvocellular pathway (P -pathway) receives input from Pp 
ganglion cells and first projects to layers 3-6 of the LGN. Cells in the parvocellular 
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Figure 2.4. Receptive field layout of a retinal ganglion cell with (a) an “ON” center 
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periphery inhibitory. From Schiffinan [1996], 
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Figure 2.5. Images were obtained from Dr. Van Essen’s laboratory home page, Washington 
University, St. Louis, Mo. 
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geniculate layers project to layer 4 C 3 in the primary visual cortex From layer 4C P they, 
in turn, project to layers 2 and 3. From these two layers, information is sent to visual area 
2, layer 4, and then to the inferior temporal cortex (Tovee, 1996). The inferior temporal 
cortex is thought to be concerned about the ‘what’ of an object. In summary, a crude but 
simple classification of the P and M pathways can be characterized as the ‘what and 
where’ of objects that an observer perceives. (Livingstone and Hubei, 1988) 

Scientists have learned about the P and M pathways and their contributions to vision 
through studies that observed one pathway after the other pathway had been made inoperative 
by lesioning it. In 1990, Schiller and Logothetis created lesions in the P or M pathway of 
monkeys and then conducted various tests, including color perception, texture perception, 
acuity, pattern perception, flicker perception, and contrast perception. The results of these 
tests indicate that particular functions do tend to be associated with a specific pathway. In 
general, the P pathway is associated with color perception, texture perception, pattern 
perception, acuity, and contrast perception, whereas the M pathway is associated with flicker 
and motion perception. While only parvocellular lesions had an effect on the monkey’s ability 
to discriminate between subtle color differences, these same lesions did not affect the monkey’s 
ability to detect a single large target whose color differed from its background, even when the 
target was equiluminous with its background. This implies that the M pathway is capable of 
conducting some gross color information and can do so at isoluminance. While lesions in the 
M pathway resulted in a deficit of flicker perception, this was universally true for high temporal 
frequencies only. For low temporal frequencies, lesions in the M pathway had no effect, thus 
demonstrating that the P pathway is capable of transmitting low temporal frequency 
information. (Sekuler and Blake, 1994) 

The P pathway also provides information to bloblike regions of the visual cortex. 
Unlike most ganglion cells, the cells in these bloblike regions are not at all concerned about 
orientation. These cells called blobs exhibit color opponency in a manner similar to that of the 
ON and OFF center cells. However, these cells turn ON and OFF in response to a specific 
chromatic illumination instead of an overall illuminance (Sekuler and Blake, 1994). Having this 
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basic physiological background, we can now focus on possible physiological explanations for a 
phenomenon known as the “oblique effect.” 

It is well documented that horizontal and vertical lines are easier to see than lines at 
oblique angles, a phenomenon known as the “oblique effect” (Campbell, Kulikowski and 
Levinson., 1966; Appelle, 1972). This phenomenon has been observed in a variety of visual 
tasks. Essock (1980) divided the oblique effect into two classes. Oblique effects arising from 
basic visual functions such as detection, acuity, and other measures of sensitivity are termed 
Class I oblique effects. The Class I oblique effect is not caused by a bias in the optics of the 
eye (Campbell and Kulikowski, 1966), but is thought to result from the orientation bias of the 
P-cells located in the visual cortex (Lennie, 1974). Several studies have shown that an oblique 
effect is most observable when a stimulus with a high spatial frequency and a low temporal 
frequency is presented to the fovea (Maffei and Fiorentini, 1973; Berkley, Kitterle and 
Watkins, 1975; Camissa, Blake andLema, 1977). 

Class II oblique effects arise from tasks that require subsequent processing of stimulus 
information. For example, a task measuring the detection threshold of a stimulus oriented 
either obliquely or non-obliquely would result in a Class I oblique effect. When an observer 
must press one of two buttons indicating whether two simultaneously-presented stimuli of 
various orientations are the same, or whether they are different, the result would be a Class II 
oblique effect. This Class II oblique effect would be the result of encoding or further 
processing of stimulus information required for task completion. An important distinction, 
however, is that the Class I oblique effect discussed above results from achromatic stimuli. For 
chromatic stimuli, the results are not clear and leave uncertainty as to whether an oblique effect 
exists (Kelly, 1975b; Murasagi and Cavanagh, 1989). This thesis explores a Class I oblique 
effect and will be specific as whether this oblique effect is a result of achromatic or chromatic 
stimuli. (Essock, 1980) 

There are various hypotheses explaining why the oblique effect exists. One suggests 
that the world we live in, especially urban areas, contains more stimuli oriented horizontally or 
vertically as opposed to obliquely; thus, visual experience plays a role in determining sensitivity 
(Annis and Frost, 1973). However, the oblique effect has been demonstrated with infants as 
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young as six weeks old (Leehy, Moskowitz-Cook, Brill and Held, 1975). It seems unlikely 
that the visual experience of infants would be sufficient to account for the oblique effect. To 
further confuse the issue, it has been shown that with extensive practice detecting diagonal 
lines, observers may improve their oblique sensitivity until it is equal to their sensitivity for 
detecting horizontal and vertical lines (Mayers, 1983; Krebs, 1992). A more widely accepted 
hypothesis suggests that more cells are tuned to vertical and horizontal stimuli than to oblique 
stimuli. 

Similarly disputed is whether the oblique effect is a retinal or a neural phenomenon. 
Campbell and Kulikowski (1966) and Mitchel and Muir (1967) used lasers to bypass the optics 
of the eye by projecting stimuli directly onto the retina. The oblique effect was obtained in 
both studies. Other studies involving head tilt (Rock and Heimer, 1957; Attneave and Reid, 
1968) further investigated whether the oblique effect was a retinal phenomenon. When 
subjects viewed tilted stimuli with their heads tilted the same amount as the stimuli, the stimuli 
were retinally upright, but phenomenally oblique. (“Phenomenal!/’ describes the stimuli’s 
orientation in the visual frame of reference of the subject [Lasaga and Gamer, 1983], The 
phenomenal frame of reference in this case would be gravitational.) Other, similar studies have 
adopted an arbitrary frame of reference (Rock and Heimer, 1957). In a study by Attneave and 
Reid (1968), subjects were told to think of the top of their heads as vertical, regardless of 
whether or not they were tilted. For head-tilt experiments in which subjects had their heads 
tilted 45 degrees, stimuli that were horizontal/vertical were retinally oblique, but were 
phenomenally upright. Unless otherwise instructed, subjects tended to adopt a phenomenal 
frame of reference rather than a retinal frame of reference. Therefore, when a subject’s head 
was tilted 45 degrees, an oblique effect was obtained for oblique stimuli, even though these 
stimuli were not retinally oblique (Attneave and Olsen, 1967). However, when subjects were 
told to think of the top of their heads as vertical, regardless of their 45 degree head tilt, they 
displayed an oblique effect for stimuli that were gravitationally upright, but retinally oblique 
(Attneave and Reid, 1968). In light of these studies, the oblique effect is highly unlikely to be a 
retinal phenomenon. 



12 



Another possible cause of the oblique effect is a neural phenomenon. This hypothesis, 
and that the origin of this neural phenomenon arises in the P pathway, seem credible even 
without the additional evidence provided by Rabin’s (1992 and 1994) work with Visual 
Evoked Potentials (VEPs). The P pathway is responsible for acuity information and, thus, 
spatial information and is capable of processing low frequency temporal information The 
Class I oblique effect has been observed primarily at high spatial frequencies and low temporal 
frequencies. 

While the origin of the Class I achromatic oblique effect has been disputed, the dispute 
about a Class I chromatic oblique effect is not over its origin, but over its very existence. One 
of the earliest articles discussing the Class I oblique effect and chromaticity is a 1975 study by 
D. H. Kelly. 

The stimuli Kelly used in the experiment were striped luminous-contrast gratings 
flickering sinusoidally. A grating is a pattern of adjacent light and dark bars or stripes, and a 
sinusoidal grating, shown in Figure 2.6, has a gradual transition from light areas to dark areas 
with no sharp edges (Schiffman, 1996). The gratings were presented horizontally, vertically. 




Figure 2.6. Sinusoidal gratings. Wide stripes correspond to low spatial frequencies. As 
the spatial frequency increases the stripes become thinner. 

and at angles of 45 and 135 degrees (obliquely). A plot of the mean threshold at each 
orientation for various temporal frequencies ranging from approximately 1-40 hertz is shown in 
Figure 2.7. As the figure shows, the thresholds for the oblique presentations are consistently 
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lower than those for the non-oblique presentations for all temporal frequencies. This result is 
consistent with previous results. (Kelly, 1975) 




Figure 2.7. Achromatic sine-wave flicker sensitivity curves. 10’ = 3 cpd and 22’ = 1.36 
cpd. From Kelly [1975], 



Although this result was consistent, Kelly wanted to test the hypothesis that luminous 
contrast was a necessary condition for the oblique effect. He repeated the experiment, but 
changed the stimulus to a red-green equiluminous grating. The thresholds obtained, presented 
in Figure 2.8, showed that for the 10-minute stripes, an oblique effect was present for temporal 
frequencies under approximately 10 hertz. However, for the 22-minute stripes, no oblique 
effect was observed. This is due to the luminous grating sensitivity decreasing with decreasing 
spatial frequency, whereas chromatic grating sensitivity remains constant (Kelly, 1975). Kelly 
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concluded that the oblique effect for the 10-minute stripes was probably a hybrid response 
resulting from a spurious luminous component. This spurious luminous component was likely 




Figure 2.8. Chromatic sine-wave flicker sensitivity curves. From Kelly [1975]. 

a result of the stimuli not being isoluminant for the observer. Obtaining exact 
isoluminance is not an easy task. (Kelly, 1975) 

Kelly provided a study that could be used directly in the dispute over a possible Class I 
chromatic oblique effect. Other studies were used indirectly and provided better tools or 
research methods. One such study by Mullen (1985) had important implications for future 
vision research in this area. 
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Mullen’s article, “The Contrast Sensitivity of Human Colour Vision to Red-Green and 
Blue- Yellow Chromatic Gratings,” described an innovation with her experimental design which 
led to measurements without any type of chromatic aberration. By using a large field size, she 
was able to measure thresholds for low spatial frequencies without the reduction in luminous 
sensitivity shown to occur with spatial frequencies below approximately four cycles/deg. 
Instead of using only one chromatic intensity value for all specific spatial frequencies, as had 
often been done before, Mullen used a number of selected points. This provided more 
accuracy. These factors combined to give the chromatic contrast sensitivity function (CSF) 
obtained for red-green gratings a much different look than previously thought. The CSF still 
had the same basic shape, but the cutoff for high frequencies occurred much earlier at 
approximately 10-12 cycles/deg. (Mullen, 1985) 

The CSF is a method used to describe the visual system’s sensitivity to sinusoidal 
waveforms. Contrast, as defined in a CSF, is a relative measure that is computed rather 
than measured. Contrast, the difference between stimuli elements, is formally defined as 
the amplitude of a waveform relative to its mean. Therefore, at a mean luminance level of 
.5 cd/m 2 , a sinusoidal grating with a contrast of 50 percent would have a trough of .25 and 
a peak of .75 cd/m 2 . This same waveform at a mean luminance level of 500 cd/m 2 would 
still have a contrast of 50 percent if its peak were at 750 and its trough were at 250 cd/m 2 
(Schiffman, 1996). The use of sensitivity (1/threshold contrast) in CSF is similar to 
everyday usage; therefore, a low detection threshold is equivalent to high sensitivity. 
(Schiffman, 1996) 

A CSF for spatial frequency is shown in Figure 2.9. Peak sensitivity is found at 
approximately three cycles per degree (cpd), with approximately 50 cpd being the cutoff 
for high frequency acuity. 

Prior to Mullen’s work, studies with red-green stimuli used frequencies more than 
20 cpd and suggested resolution above 25 cpd (Mullen, 1985). Mullen’s work provided a 
template for further research. It is not by coincidence that Murasagi and Cavanagh’s 1989 
article dealing with the chromatic oblique effect, and using red-green stimuli, chose spatial 
frequencies under the 10-12 cpd cutoff proposed by Mullen. 
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This article by Murasagi and Cavanagh further explored earlier work by Kelly 
regarding the oblique effect for luminous, as well as chromatic, stimuli. Kelly had 
postulated the absence of an oblique effect for chromatic channels if the opponent-color 




Figure 2.9. Contrast Sensitivity Function. From [Schiffrnan, 1996], 

pathways for humans, like those in monkeys, were not orientation selective. However, 
research published the same year as Kelly’s article (Poggio, Baker, Mansfield, Sillito and 
Grigg, 1975), as well as additional research a few years later (Michael, 1978), revealed 
that monkeys might possess orientation selectivity in their chromatic channel. The 
possibility of the chromatic channel analyzing orientation independently of the luminous 
channel led the authors to design an experiment to test this possibility in humans by 
determining if an oblique effect obtained with chromatic stimuli differed from that 
obtained with strictly achromatic stimuli. (Murasagi and Cavanagh, 1986) 

To test this possibility, the researchers used a constant temporal frequency of 2 Hz 
and spatial frequencies of 2, 4 and 8 cpd were used. The stimuli were sinusoidal gratings. 
The gratings were presented at oblique (45 and 135 degrees) and non-oblique (0 and 90 
degrees) angles. Axial chromatic aberration was taken into account by having the subjects 
view the stimuli through an achromatizing lens. A revised ascending method of limits was 
used to determine thresholds for both luminance and chromatic stimuli Since the 
production of an isoluminant stimulus is a non-trivial matter, with isoluminance varying 
slightly from subject to subject, Murasagi and Cavanagh used stimuli in five different areas 
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in the neighborhood of equiluminace. The maximum threshold for the spatiotemporal 
region they were investigating was assumed to occur at equiluminance. They made this 
assumption because a chromatic grating with no luminous component should be the 
hardest to detect, as detection is by color alone rather than by color and luminance. 

A significant main effect of orientation was present for all observers, as was a 
significant three-way interaction between grating types. Spatial frequency and orientation 
effects were present for three of the four observers. This showed that, for three of the 
four observers, the effects of the four orientations at certain spatial frequencies were 
different for achromatic and chromatic stimuli. 

Like Kelly, Murasagi and Cavanagh have possible problems with spurious 
luminous components. By taking five measurements in the neighborhood of 
equiluminance, the contribution of these components has probably been reduced. 
However, if the actual isoluminance point were in a region between the areas they picked, 
then a luminous component would be present in their stimuli. Additionally, while spatial 
frequencies of 2, 4 and 8 cpd were used, a chromatic oblique effect was present only at 8 
cpd for three of the four observers. At 8 cpd, chromatic aberration is a factor. The 
researchers used an achromatizing lens to account for axial chromatic aberration, but they 
did not take into account lateral chromatic aberration. Post hoc analysis minimized the 
possibility of this being a factor. As the authors themselves state, any slight misalignment 
between a subject and the achromatizing lens would result in a spurious luminance 
component. (Murasagi and Cavanagh, 1986) 

During the same time frame as the Murasagi and Cavanagh work, Bradley, 
Switkes and De Valois, were also exploring Kelly’s earlier work. The authors designed an 
experiment to compare the visual processing of chromatic and luminance information. 
The prolonged viewing or adaptation of a sinusoidal grating desensitizes the observer to 
similar gratings, especially when the similarity is in orientation or spatial frequency. 
However, this desensitivity is not present for gratings with orientations differing by 
approximately 45 degrees or spatial frequencies differing by 1.5 octaves. Thus, this 
adaptation has been termed selective. 
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The effects of selective adaptation have been used as psychophysical evidence for 
the presence of spatial frequency-selective and orientation-selective neurons in the human 
visual system. However, the behavior of cells displaying selective adaptation for spatial 
frequency when measured psychophysically has not always been consistent with their 
physiology. Thorell, De Valois, and Albrecht (1984) observed neurons that displayed 
different spatial frequency tuning depending on whether the stimulus contained luminance 
or color. They observed low-pass tuning for chromatic gratings, but band-pass tuning for 
luminance gratings. Additional studies by Livingstone and Hubble (1984) and Lennie, 
Sclar and Krauskopf (1985) found that cells in the visual cortex that responded to 
isoluminant color contrast did not display selective adaptation for orientation or spatial 
frequency. (Bradley, Switkes and De Valois, 1988) 

Bradley et al. (1988) set out to explore this inconsistency with spatial frequency 
adaptation and orientation for chromatic gratings. The zero contrast condition for all 
gratings was a uniform yellow field with a chromaticity that was adjusted for each 
observer’s differing sensitivity to red and green phosphor emissions. This varying of the 
zero contrast condition enabled presentation of the red-green sinusoidal gratings at each 
observer’s isoluminance axis. Both isochromatic and isoluminant gratings were presented 
and were viewed through an achromatizing lens. To overcome problems associated with 
making repeated measurements of a decaying effect, the researchers used a long initial 
adaptation period followed by alternation of a brief stimulus presentation with a brief 
adaptation period. The stimuli used for adapting was a 2 cpd grating, run separately for 
each of the four possible conditions of horizontal or vertical and luminance or chromatic. 
For a spatial frequency of 2 cpd, thresholds for luminance gratings were similar in both 
pre- and post-adaptation trials for oblique angles and showed the desensitivity expected at 
horizontal and vertical angles. However, while the pre- adaptation data for chromatic 
gratings at this frequency did not show an oblique effect, an oblique effect was evident in 
the post-adaptation data. (Bradley, et al., 1988) 

A similar experiment varying spatial frequency while keeping orientation constant 
confirmed that, for varying spatial frequencies, a specific spatial frequency adaptation 
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effect can be observed for sinusoidal luminous gratings. This experiment was repeated 
with chromatic gratings, and a specific spatial frequency adaptation effect was also 
observed with sinusoidal chromatic gratings. The results of this study along with previous 
psychophysical studies demonstrating the parallels between the data for luminance and 
color (De Valois and Switkes, 1983; Switkes and De Valois, 1983; and Ware and 
Mitchell, 1974) suggest that for the beginning stages of human vision, color and luminance 
are processed in a similar manner. (Bradley, et al., 1988) 

In addition to the numerous psychophysical studies on the oblique effect, other 
studies have been conducted electrophysiologically. When studying primates or other 
animals, collecting data is often not possible through psychophysical means. Although an 
animal may not be able to verbalize or react, a response may still be obtained 
electrophysiologically. By electrically stimulating an individual cell, it is possible to 
monitor the cell output. VEPs provide an additional method of studying the role of 
chromatic patterns in perception (Rabin, 1994). In Rabin’s 1992 paper “VEP’s in Three- 
Dimensional Color Space,” a Class I oblique effect at isoluminance or a chromatic oblique 
effect was shown at the spatial frequency of 1 cpd. Psychophysically, the Class I oblique 
effect for luminance or chromaticity is typically not obtained at low spatial frequencies. 

The Class I oblique effect for achromatic stimuli has been obtained under a number 
of different conditions. It has been demonstrated psychophysically (Campbell and 
Kulikowski, 1966; Camisa, et al., 1977) and electrophysiologically (Maffei and Campbell, 
1970; Rabin, Switkes, Crognale, Schneck and Adams, 1994). Electrophysiologically the 
oblique effect is evident by comparing the output of microelectrodes for oblique and non- 
oblique stimuli. These microelectrodes monitor the VEPs of cells as they are exposed to 
different stimuli. Psychophysically, the Class I oblique effect is evident by comparing the 
responses of subjects for oblique and non-oblique stimuli. A Class I chromatic oblique 
effect could be measured the same way. However, the Class I oblique effect for chromatic 
stimuli has not been obtained under various conditions, and whether such an oblique effect 
actually exists is a matter of debate. To participate in this debate, it is necessary to 
understand how the information that the eye collects is processed. One explanation is that 
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the eye uses a process similar to Fourier analysis so named after the nineteenth century 
mathematician responsible for this analysis. 

Jean Baptiste Fourier studied how heat flows through an object when it is heated 
up and found that heat behaved in waves. He modeled these waves using complex 
equations, and discovered that they consist of periodic waveforms. Fourier found that any 
quantity that changed in a complex manner over time could be converted into a series of 
simple sinusoidal functions. Each sinusoidal pattern could be defined by its period, 
frequency, and angular velocity. This process is now known as Fourier analysis. (Who is 
Fourier?, 1995) 

Fourier analysis can be used to analyze a natural scene by decomposing it into a 
sum of a series of sinusoidal components, each having a different spatial frequency, 
amplitude, and orientation. Vision scientists believe that the human visual system uses a 
process similar to Fourier analysis to process visual imagery. The human eye receives 
different intensities of light reflected from an object. These light intensities pass through 
the cornea and filter down into the photoreceptors. The photoreceptors then send an 
electrical signal to the brain, where these neural responses are categorized into specific 
spatial channels. Psychologists believe that this sensory input is transformed into a neural 
response, which is then categorized into a perceptual experience. If the visual system 
passes the image, and this image corresponds to a perceptual experience, then the observer 
can recognize the object. However, if your cornea is degraded— e.g., a cataract— high 
spatial frequency sinusoidal waves will not pass through the lens and be sent to the brain. 
A degraded signal such as this or a lack of a similar perceptual experience may result in 
failure to recognize the object. The absence of high spatial frequencies will cause the 
image to appear blurry. In some cases, the amplitude of those missing high spatial 
frequencies can be increased so that these signals can be sent to the brain. 

Scene (Figure 2.10) can be broken down into many different visual components by 
Fourier analysis or other tools. These components or parameters include color, 
orientation, and spatial qualities and include the scene as a whole, as well as for the 
individual objects that comprise the scene. By manipulating the parameters of an object 
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(e g., armored personnel carrier [APC]) in a scene, it is possible to camouflage this object. 
This may be done by changing the object’s color to match that of its background through 




Figure 2. 10. Armored Personnel Carrier 



temporary means such as netting, or by more permanent means such as paint. A Fourier 
analysis shows the spatial composition of the scene. The low spatial frequencies are 
located in the center of Figure 2.11, and the high spatial frequencies are found in the 
comer regions of the figure. The high spatial regions result from the edges of the APC. 
These high spatial frequencies contrast with the low spatial frequencies found elsewhere in 
the scene. Netting would reduce these high spatial frequencies and would also lessen the 
edge effect evident in Figure 2. 12, thereby enhancing the APC’s camouflage. 

We have looked at a scene’s color, orientation, and spatial information and how 
these parameters can be manipulated to achieve better camouflage. The parameters 
manipulated in the experiment in this thesis are spatial and orientation. Since the Class I 
oblique effect has been primarily observed at low temporal and high spatial frequencies, a 
temporal frequency of 2 Hz and spatial frequencies of 1, 3 and 7 cpd were chosen. The 
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11. Fast Fourier Transformation on the APC. 
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Figure 2. 12. High-pass filter of the APC originally shown in Figure 2.10. 



red-green spatial CSF begins to decrease at frequencies greater than 1 cpd; therefore the 
frequencies of 3 and 7 were chosen knowingly, trading off sensitivity for the advantages of 
a higher spatial frequency where there would be a higher likelihood of observing an 
oblique effect. Frequencies higher than 7 cpd were not chosen due to increasing effects of 
chromatic aberration. 

A Class I chromatic oblique effect was expected to be observed at spatial 
frequencies of 3 and 7 cpd. Psychophysically, the Class I oblique effect has not been 
readily observed at spatial frequencies as low as 1 cpd and accordingly a Class I chromatic 
oblique effect was not expected to be observed at this spatial frequency. 
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III. METHODS 



A. SUBJECTS 

The experiment was conducted concurrently at the Naval Postgraduate School 
(NPS) and the University of Louisville, Kentucky (UL). Four subjects, 2 NPS and 2 UL, 
volunteered for this experiment. All subjects had normal (20/20), or corrected to normal, 
acuity and color vision. Color vision was verified with pseudo-isochromatic plates. Two 
of the four subjects (1 NPS and 1UL) were naive as to the purpose of the experiment. 
The other two subjects, the author and the remaining UL subject, were experienced 
psychophysical observers. All subjects signed an informed consent and were briefed on 
the ethical conduct for subject participation specified in the Protection of Human Subjects, 
SECNAV Instruction 3900. 39B. Subjects were screened for uncorrected astigmatic 
errors by determining spatial resolution limits for 0°, 45 °, 90 °, and 135 °. 

B. APPARATUS 

Stimuli were presented by a VisionWorks computer graphics system (Vision 
Research Graphics, Inc.) on an IDEK MF-8521 high resolution color monitor (21" X 20" 
of viewable area) equipped with an non-glare, anti-reflect, P-22 phosphor. The monitor 
had a resolution of 800 by 600 pixels (x=75.02 and y=74.92 pixels/degree), 98.9 Hz 
frame-rate, mean chromaticity of r = 0.334 , g = 0.336, b = 0.300 (1931 CIE), and a 
maximum luminance of 100 cd/m2. Refer to Table 3.1 for the chromaticity and 
luminance coordinates for each phosphor. The University of Louisville’s apparatus and 
procedure were identical to the Naval Postgraduate School’s, except that the stimuli were 
displayed on a 17” Nanao Flexscan F2.21 color monitor. Subjects viewed the monitor 
from 1.5 meters and were positioned by an adjustable chinrest. A small floor lamp (2 6 
cd/m2) was positioned behind the monitor to reduce screen glare. 
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cm 





X 


V 


z 


Luminance (cd/m 2 ) 


Red phosphor 


.617 


.345 


.038 


24.0 


Green phosphor 


.334 


.581 


.085 


88.7 


Blue phosphor 


.162 


.081 


.757 


12.7 



Table 3.1. Chromaticity and luminance of monitor 
C. STIMULI 

Sinusoidal gratings were presented within a spatially windowed circular test field that 
subtended 7.59° of visual angle. The Gaussian window was truncated at ±1 standard 
deviations for both x and y directions. The test patterns were one-dimensional spatio-temporal 
sinusoids of varying orientation (principal and oblique), spatial frequency (1.0, 3.0, and 7.0 
cycles/degree), and color contrast. Test patterns for each subject consisted of two orientations, 
principal (0° and 90°) and oblique (45° and 135°). For each subject, maximum sensitivity for 
each orientation within the principal and oblique grouping was chosen. All sinusoids were 
raised cosines temporally modulated at 2.0 Hz. The sinusoid pattern was presented in a 
1500 msec interval with contrast ramped on and off according to a linear window. 
(Contrast peaked at 202 msec and fell at 1304 msec). 

Color contrast was computed by different ratios of percent red and green (Sellers 
et al.,1986). The monitor was controlled by a Cambridge Research Systems VSG 2/4 
video board that was linearized to 10 bits of resolution per gun. The outputs of each gun 
were linearized by means of stored look-up table file. Sixteen different sinusoidal red- 
green color mixtures were generated by changing the red phosphor only, green phosphor 
only, or by changing the red and green guns in fixed proportions. Color contrast was 
defined according to the (Michelson) formula shown in Equation 3.1. Blue gun was held 
constant in all quandrants. Red and green gun values were used in the determiniation of 
red and green contrast as shown in Equations 3.2 and 3.3. 
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contrast = 


(peak - trough) / (peak + trough) 


3.1 


red contrast = 


(red gun value - 50) / 50 


3.2 


green contrast = 


(green gun value - 50) / 50 


3.3 


D. PROCEDURE 


Thresholds were 


determined by a two-altemative forced 


choice adaptive 



psychometric procedure, QUEST (Watson and Pelli, 1983). Threshold was defined at 75 
percent correct. A total of 480 trials, 30 trials per condition, were randomly presented 
within each session. A session (~ 45 minutes) consisted of one sinusoidal condition with 
16 different red-green color mixtures. A subject' had to complete six sessions to 
contribute one threshold point for all conditions. 

At the beginning of each session, subjects dark-adapted for approximately five 
minutes before initiating the first experimental trial. Three of the four subjects were tested 
monocularly, while the fourth subject (UL) was tested binocularly. At the beginning of 
each trial, the subject was instructed to focus on a fixation cross (.19° by .13°) located in 
the center of the screen. The subject initiated the first trial with a keyboard response, the 
fixation cross extinguished followed by presentation of the first interval, 121 msec ISI, and 
then presentation of the second interval. The subject’s task was to detect which interval 
contained the sinusoidal grating. The next trial followed 250 msec after the subject’s 
keyboard response. 

Color contrast thresholds were determined from 16 different color-mixture ratios. 
The sixteen different ratios could be divided into four different percent red and green 
quadrants. Quadrant one started with 100 percent green and 0 percent red, quadrant two 
started at 100 percent red and 0 percent green, quadrant three started at -100 percent 
green and 0 percent red, and quadrant four started at -100 perecent green and 0 percent 
red. The red-green ratios within each quadrant were 0.0, 0.5, 1.0, and 2.0. The 
thresholds from these red-green ratios will form an ellipse with the half-length of the axis 
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measuring color discrimination and the half-width of the axis measuring brightness 
discrimination (Sellers, 1986). 
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IV. RESULTS 



Many researchers have carried out experiments to determine visual sensitivity to color 
differences. One way of determining these differing sensitivities is through differing values of 
International Commission on Illuminance (ICI or, more commonly, CIE, for the French translation) 
primaries (Kaiser and Boynton, 1996). These primaries allow all colors to be specified in terms of 
three numbers representing the red, green and blue primaries. The color-matching type experiment 
is set up to test if an observer can discriminate between a chosen color and another color similar to 
this color. Color-matching experiments have looked at the standard deviations of color-matchings 
for representative colors throughout the color spectrum. These standard deviations are directly 
related to the corresponding just-noticeable difference of colors. (Brown and MacAdam, 1949) 

Numerous surveys of differential thresholds have been carried out, but W. D. Wright (1941) 
and D.L. MacAdam (1942) completed two of the more extensive surveys. “MacAdam plotted the 
results of the survey on the chromaticity diagram in terms of the standard deviation of color-matching 
in several directions for selected colors.”(Brown and MacAdam, 1949) The figures resulting from 
this survey formed closed curves on the diagram and the closed curves were elliptical in shape. 
(Brown and MacAdam, 1949) 

In a later paper, Silberstein and MacAdam discussed that the errors of these closed curves 
were Normally distributed. They deduced that the curves should be ellipses, as they appeared to be. 
They further expostulated that if the variations were not confined to chromaticity, the closed curves 
would form ellipsoids rather than ellipses. Using the assumption that the probability of making a 
match that falls within a specific region of color space, near the target color, was not changed by any 
change of primaries, Silberstein proved the standard deviation figures to be ellipsoids. (Brown and 
MacAdam, 1949) 

The fact that there was a theoretical explanation for MacAdam ’s ellipses and not just an 
empirical observation was interesting. However, the discrimination ellipsoids of Brown and 
MacAdam were obtained for a bipartite field only. Noorlander, Heuts and Koenderink (1979 and 
1980) and Noorlander and Koenderink (1983) furthered the work of Brown and MacAdam by 
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extending their methodology from a simple bipartite field to more complicated stimuli by varying 
temporal and spatial frequencies. When three primaries are used discrimination ellipsoids are 
obtained. However, if just two primaries are used, such as red and green, than a cross-section of a 
discrimination ellipse is obtained. Neurlander, Heuts and Koenderink (1980) did this by obtaining 
discrimination ellipses for a number of different spatial and temporal frequencies. The lengths of the 
major and minor axes of these ellipses, as well as their orientation, are highly dependent on both the 
spatial and temporal frequencies. 

The finding that the cross-section of a discrimination ellipsoid is an ellipse was also used by 
Sellers et al. (1986) in a study of congenital and acquired color defects. However, in the Sellers et 
al. paper, the axes of their graphs were not primaries as in Brown and MacAdam’s or Neurlander , 
Heuts and Koenderink’ s, but were percent red contrast and percent green contrast. Even without 
knowledge of discrimination ellipsoids, the fact that Sellers et al.’s data form an ellipse can be 
explained by examining the following model. The central dashed line in Figure 4.1 is the 
equiluminance axis. Assume there are two luminance processes for detecting the brightening and 
darkening of a spot. The thresholds of these processes are displayed in Figure 4. 1 with dashed lines 
and are labeled “BRIGHT” and “DARK.” Similarly, the processes for detecting color are labeled 
“RED 1” and “GREEN 1.” (“RED2" and “GREEN2" refer to two different thresholds.) As a first 
approximation, the visual threshold will be determined by whichever process has the lowest threshold. 
Therefore, this will be the parallelogram bounded by the four lines “BRIGHT,” “DARK,” “RED1," 
AND “GREEN1.” A phenomenon known as probability summation accounts for the rounding of the 
comers of the parallelogram, and an ellipse is formed. Probability summation occurs near the comers 
of the parallelogram and can be thought of as a sum of two processes, e.g. “BRIGHT” and “RED1 ." 
For example, if the probability of either of the processes detecting a stimulus is .5, then if both 
processes are independent, the probability of either of them (or both) detecting a stimuli is .75. Thus, 
a contour connecting points where the probability of detection is .5 will exclude the comers of the 
parallelogram (Graham, 1989). (Sellers, 1986) 

A ratio can be determined by dividing the length of the major axis by the length of the minor 
axis. For ellipses with a ratio greater than four, the major axis nearly coincides with the 
equiluminance axis. The angle of discrepancy between the equiluminous axis and the major axis 
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Red-green mixture thresholds 



GREEN 1 



8R.ICH T 



:pARK: 



RED. CONTRAST, % 



Figure 4.1. Detection model. From Sellers et al.[1986]. 



varies approximately as the inverse square of the length/width ratio. Therefore, elongated ellipses 
have a major axis that nearly coincides with their equiluminous axis. (Sellers et al., 1986) 

The paper by Sellers et al. (1986) is extremely important to this thesis in that its methodology 
provided a foundation for the methodology used in this thesis. Using percent red contrast as the x- 
value and percent green contrast as the y-value, thresholds were determined for 16 different rays. 
For each ray, the proportion of percent red contrast to percent green contrast is constant along the 
ray. When plotted, the thresholds form an ellipse where the half-length of the major axis is a useful 
measure of color discrimination, and the half- width is a useful measure of brightness discrimination. 
Sellers et al. were interested in length and orientation, since they used these values for classification 
of color deficient subjects. Major axis length and orientation are important in this thesis. However, 
the crucial fact that Sellers et al.’s methodology resulted in data that theoretically forms ellipses is the 
crucial item. Data points from rays in the vicinity of isoluminance will have high leverage since they 
will be close to the end of the major axis, but they do not need to be at isoluminance. 
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The subjects nan in this experiment did not have any color defects. Thus, collecting data on 
these subjects for classification purposes was not an exceptionally interesting endeavor. However, 
the possibility of using this methodology to explore a chromatic oblique effect was interesting. If 
oblique and non-oblique sensitivities are the same, and if oblique and non-oblique information is 
processed in an identical manner, then the ellipse obtained from a subject responding to non-oblique 
(horizontal or vertical) chromatic stimuli and the ellipse obtained from the same subject at the same 
temporal and spatial parameters, but with oblique chromatic stimuli, should be identical. A 
“spurious” luminous component is not a problem. Since the data points theoretically form an ellipse, 
the requirement for a point exactly at isoluminance no longer holds. 

The elliptical nature of the data has been used to fit ellipses to the data by the method of 
maximum likelihood. A well-known result from linear regression informs us that the method of 
maximum likelihood is identical to the method of least squares in this case (Larsen and Marx, 1986). 
The programs used here to fit ellipses minimize the sum of the squared error. Assistant Professor 
Professor Samuel Buttrey of the Naval Postgraduate School in Monterey, CA. created these 
programs, their sub-programs and other programs of use. He created these programs, which are 
found in Appendix A, using the statistical package S-Plus. 



The following terminology will be used to describe ellipses and their parameters (Figure 4.2). 



terminology 


definition 


a 


half-length of the major axis of an ellipse 


b 


half-length of the minor axis of an ellipse 


0 


angle as measured from the x-axis to the major axis of an ellipse 


*/ 


x-coordinate for data point i. X-axis is percent red contrast 


y, 


y-coordinate for data point i. Y-axis is percent green contrast 




x-coordinate for the center of an ellipse 


y e 


y-coordinate for the center of an ellipse 


r, 


distance from (x c , y c ) to a point on the ellipse along the ray 


e, 


error term 


U 


polar angle of r i 
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Five-Parameter Ellipse 




Figure 4.2. Five parameter ellipse. Created by Professor Samuel Buttrey. 

The parameters for the true ellipse are unknown, but they may be estimated from the data. Three 
models were used to represent the underlying ellipse, and ellipses will be classified according to the 
method in which they were modeled. The polar angle relative to the coordinate t t is fixed at the start 
of the experiment and is determined by tan(r^) =y/x i . The sixteen values for t, in degrees are (0, ±30, 
±45, ±60, ±90, ±120, ±135, ±150,180). Each of theses ellipses possesses an equiluminous axis along 
which, by definition, luminance is constant. The exact determination of this equiluminous axis is 
difficult, but for ellipses with an a/b ratio greater than four, this equiluminous axis is closely 
approximated by the major axis of the ellipse (Sellers, et al., 1986). A ray that coincides with the 
equiluminous axis will vary in color along r t , but will have a constant luminosity. Any ray that is not 
aligned with the equiluminous axis will have a constant red to green percent contrast ratio along the 
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ray, but luminosity will not be constant along the ray. Both x and y values can be positive or 
negative. 

A. CLASS I-TYPE ELLIPSES 

The Class I-type ellipse has five estimated parameters. These parameters are a, b, 0, x 0> and 
y Q . The predicted r ;Or A, is a function of the estimated parameters and r. =f[a,b,d, x c y c t\. The 
actual model used is r i = r j + €j. Here ^ are independent identically distributed (iid) and N(0,o 2 ). 

The function f{a,b,Q, x c y c 0 is complicated and can be found in Appendix A (ell.pred). The 
objective function is the sum of the squared differences between the observed r t and the predicted 
r . . Data were collected from both oblique and non-oblique stimuli; the class I-type ellipse is an 
ellipse that is obtained by fitting an ellipse to all of the data for a specific spatial frequency. For 
example, Subject One completed five runs at each condition. Each run results in the calculation of 
a threshold along each ray; thus, 16 thresholds were determined for this subject on five different 
sessions. This resulted in 80 data points for the non-oblique condition and 80 data points for the 
oblique condition for each of the three spatial frequencies used. Class I-type ellipses are fitted to the 
combined data of a subject at a specific spatial frequency. For Subject One, 160 data points were 
used, and an ellipse was fitted by the method of maximum likelihood. In the past, the ellipses 
obtained in this manner have been forced to have their center at the origin (Sellers, et al., 1986). 
However, much better fitting ellipses are obtained by allowing the center not to be pinned to the 
origin. The centers obtained for most subjects were generally close to the origin. The fact that a 
better fit was obtained by letting the ellipse be centered at coordinates other than the origin may be 
an indicator that centers may have some sort of bivariate distribution across subjects, or it may be an 
effect caused by the monitor. However, data from UL were definitely not centered on the origin and 
tended to have a center in the second quadrant. The ellipse programs in Appendix A allow the ellipse 
center to be pinned at the origin (fit.center=F) or to “float” (fit.center=T), but it did not make sense 
to pin the center to the origin when some of the actual ellipse centers obtained were definitely not at 
the origin. 
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B. CLASS-n TYPE AND CLASS-DI TYPE ELLIPSES 



The Class-II Type ellipse has six (program ellipse.!!) and the Class-IH Type ellipse has seven 
(program ellipse.III) estimated parameters. The models for the r for the three classes of ellipses are 
shown below. The model for Class-13 Type ellipses may be changed so that it is f(a, b + 6 b , 0, x e 
and yj by changing which.type from one to two in the program ellipse.!!. 



f {class I type ellipses) 


= Mb, Q,x c ,y c ) 


4.1 


r {class II type ellipses) 


= A a+ 5 a ,b,Q,x c ,y c ) 


4.2 


r {class III type ellipses) 


= Aa + ^> a ,b+6 b ,d,x c ,y c ) 


4.3 



The Class-II Type and Class-HI Type ellipses use the information of whether the data were 
from an oblique or non-oblique condition. An ellipse is then fitted to the data, but the additional 
information of whether the data were from an oblique condition or a non-oblique condition is used 
to determine a 6 a and/or a 6 b . This is done by fitting a Class-II or Class-Ill Type ellipse to the data. 
Actually, two ellipses are fit to the data, one to the non-oblique data and one to the oblique data; but 
a common ellipse center and theta are maintained for both ellipses. If the sum of the objective 
functions for the ellipse fitted to the non-oblique data and the ellipse fitted to the oblique data is 
smaller than the objective function for the ellipse fitted to all of the data, then there will be a 6 a and/or 
a 6 b . If the 6 a or6 b are small (This will be quantified in the next section.), then this orientation 
information does not significantly improve the fit of the ellipse. If the 6 a and 6 b are large, then this 
additional information significantly improves the fit of the ellipse. 



C. STATISTICAL TESTS 

A complete discussion of the statistical test used for comparing ellipses can be found on pages 
103-104 of Bates and Watts (1988). This test is an approximation, due to non-linearity, of an F-test. 
The derivation of this F statistic and how it is obtained are shown in Table 4.1 and Equation 4.4 
respectively. For this experiment, the full model consists of Class-Ill type ellipses, and the partial 
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model consists of Class-II type ellipses. The full and partial models were chosen this way because 
of interest in the significance of the difference in the length of the major axis between non-oblique and 
oblique data. 



Source 


Sum of 
Squares 


Degrees of 
Freedom 


Mean Square 


F Ratio 


Extra parameters 
Full model 


S e =S p -S f 

Sf 


v e = P f -P p 
v f = N-P f 


S e 2 = S e /V e 

Sf S f / Vf 


Se 2 / Sf 2 


Partial Model 


Sp 


N - P p 







Table 4.1. Extra sum of squares analysis. From Bates and Watts [1988], 



s / v f 



approximately 



F. 



vjN-P. 



4.4 



If the equiluminous axis is identical to the major axis of an ellipse, then if the length of the 
major axis for oblique stimuli is greater than the length of the major axis for non-oblique stimuli, a 
chromatic oblique effect has been observed. Additionally, if the length of the ellipse minor axis for 
oblique stimuli is greater than the length of the minor axis for non-oblique stimuli, a luminous oblique 
effect has been observed. For ellipses with a major axis to minor axis ratio of four or greater, the 
equiluminous axis is closely approximated by the major axis (Sellers et al., 1986). The ellipses 
obtained generally had a ratio less than four, but the major axis was still used as the equiluminous axis 
as a rough approximation. 

Four subjects took part in this experiment. Data from Subjects One and Three were collected 
at the NPS, while data from Subjects Two and Four were collected at UL. The data for these 
subjects at a spatial frequency of one cpd and both oblique and non-oblique orientations are shown 
in Figure 4.3. The measurements for a, b, theta, x c and y c are shown on the individual graphs and are 
displayed in Table 4.2, as well. 
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Sub 1 1 cpd 2 hz 



non-oblique = points/solid line oblique = pluses/dashed line 




Sub 2 1 cpd 2 hz 



non-oblique = points/solid line oblique = pluses/dashed line 




percent rad contrast 

non-cbtique e» 0.04703 b- 0.03424 cbbque » 0.04691 b= 0.0351 
thete* -1 527 x* 0.0045 y- 0.0006 p-v*Jue» 0.975 



Sub 3 1 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 




percent rad contrast 

norvodrqu* «= 0.05541 b* 0.02196 0 04845 b« 0.01964 

thete* -11.349 -0.0028 y* 0.0057 p-v»lu*« 0.091 



Sub 4 1 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 




percent rad contrast 

non-ofctrque e> 0.06444 b* 0.04976 oblique a= 0.06013 b« 0.04103 
thet»= 36.733 x* 0.0059 y= -0.004 p-v»lue* 0512 



percent rad contrast 

nan-obfiqu* e» 0.09675 b* 0.02294 obfcque e- 0.09061 b- 0.0304 
theta- -15.689 x= -0.0033 y* 0.0053 p-v»lue« 0.513 



Figure 4.3. Ellipses and values for Subjects 1-4 at 1 cpd 



Ellipse Type 


Subj 


a 


b 


theta 


center jc 


center.y 


delta .a 


delta.b 


*p-value 


non -oblique 


1 


0.04703 


0.03424 


-1.5274 


0.00447 


0.00064 








oblique 


1 


0.04691 


0.0351 


-1 .5274 


0.00447 


0.00064 


-0.0001 


0.0009 


0.975 


non-oblique 


2 


0.05541 


0.02196 


-11.349 


-0.0028 


0.00569 








oblique 


2 


0.04845 


0.01964 


-11.349 


-0.0028 


0.00569 


-0.007 


-0.0023 


0.091 


non -oblique 


3 


0.06444 


0.04976 


36.733 


0.00593 


-0.004 








oblique 


3 


0.06013 


0.04103 


36.733 


0.00593 


-0.004 


-0.0043 


-0.0087 


0.312 


non-oblique 


4 ' 


0.09675 


0.02294 


-15.689 


-0.0033 


0.00525 








oblique 


4 i 


0.09081 


0.0304 


-15.689 


-0.0033 


0.00525 


-0.0059 


0.0075 


0.513 


*Ho delta. a=0 




















Ha delta. a<>0 ! 



















Table 4.2. Ellipse values for Subjects 1-4 at 1 cpd 
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The p-values were obtained by the F-test approximation discussed in the models section; 
however, an example is shown below for further clarification. For Subject One the ellipse from the 
non-oblique data and the ellipse from the oblique data are almost identical, whereas for Subjects Two 
and Three the opposite of a chromatic oblique effect (oblique ellipses with shorter a’s than non- 
oblique ellipses) is displayed and for Subject Four, a chromatic oblique effect is shown. However, 
with an alpha of .05, none of these results are significant. 

Here is an example of how the p-values were calculated. Subject 1 completed 5 runs at all 
conditions. For each run, a total of 16 thresholds were calculated, so for a spatial frequency of 1 cpd 
a total of 80 oblique data points and 80 non-oblique data points were collected for this subject. The 
data was input to the program ellipse.III, and the ellipse center was allowed to float or not be pinned 
to the origin. From the output of this program, the objective function value is obtained. This 
objective function value is the sum of the squared error (SS m for Class-Ill Type ellipses). This value 
is subtracted from the objective function value obtained from the output of a Class-II Type ellipse 
obtained with the same data, SS n . The difference is then divided by the difference in the number of 
estimated parameters, or degrees of freedom, between the two classes of ellipses. This is the 
numerator for the equation. There is only one additional estimated parameter for Class-Ill Type 
ellipses, compared to Class-II Type ellipses, so this number is a one. Finally, the denominator is the 
value of the objective function obtained from the Class-Ill Type ellipse output divided by its degrees 
of freedom. For Class-Ill Type ellipses, seven parameters are estimated, so the 160 degrees of 
freedom for subject One decreased to 153 degrees of freedom. The resulting fraction is shown in 
equation 4.5. This fraction is referred to the F distribution with 1 and 153 degrees of freedom. Thus 
the fraction is approximately an F random variable with degrees of freedom given by 1 and 153, 

(SSn - SSmW) 45 

SS m / 153 

if the hypothesis that the ellipses differ only by 6 a , and if iid Normal errors are true. A p- value is then 
calculated through tables or statistical programs. P-values for other subjects were calculated 
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similarly, with the degrees of freedom reflecting the number of observations the subject had for that 
condition. 

The data for subjects at a spatial frequency of three cpd and both oblique and orientations are 
shown in Figure 4.4. The measurements for a, b, theta, x c and y c are shown on the individual graphs 
and are displayed in Table 4.3, as well. A chromatic oblique effect is shown for Subjects One, Two 
and Three and is significant for Subjects One and Two. 

The data for the subjects at a spatial frequency of seven cpd and both oblique and non-oblique 
orientations are shown in Figure 4.5. The measurements for a, b, theta, \ and y c are shown on the 
individual graphs, and are displayed in Table 4.4. An achromatic oblique effect is expected here and 
is evidenced by 6 b <>0 leading to larger b values for oblique ellipses compared to the b values for non- 
oblique ellipses. A peculiarity of the program that determines the 6 values is that if the 6 value is 
negative and if it is larger in magnitude than the value to be added to, then the signs of both values 
must be reversed. The hypothesis of 6 b <>0 was tested in a manner similar to 6 a <>0, and all subjects 
displayed an achromatic oblique effect. This achromatic oblique effect was significant for Subjects 
One, Two and Four. A chromatic oblique effect is shown for Subjects One and Two, but is not 
significant for either. 

The data collected from the subjects were extremely variable. This variability is not only from 
subject to subject, but also from day to day and run to run. To display some of this variability. Table 
4.5 shows p-values (uncorrected for multiple comparisons) for a subject’s run at a specific condition 
against all of the other runs at this identical condition. 

In summaiy, at one cpd neither an achromatic nor a chromatic oblique effect was shown. At 
three cpd a chromatic oblique effect was shown for three subjects and was significant for two of 
them. At seven cpd both achromatic and chromatic oblique effects were shown. All four subjects 
showed an achromatic oblique effect and this oblique effect was significant for three of them. Only 
two of the four subjects showed a chromatic oblique effect, but neither of the p-values were 
significant. 
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P#rc»nt clean contrast 



Sub 1 3 cpd 2 h z 

non-oblique = points/solid line oblique = pluses/dashed line 




Sub 3 3 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 





1 . ^ a 






i 

♦ 


• - » - • 



•0.10 -0.05 0.0 0 05 0.10 



Sub 2 3 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 




parcant rad contrast 

non-obGquo »» 0.0459 b= 0.01365 obliqua a* 0.0556 b“ 0.01 564 
that* = -16 952 x* -0.0153 y« 0.0131 p-valua- 0.005 



Sub 4 3 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 




- 0.2 - 0.1 0.0 0.1 0.2 



pore ant rad contrast 

non-obtrqua »= 0.08941 b= 0.02025 ottiqja a= 0.10976 b= 0.02385 
that*- -1 2315 x» 0 008 y» 0 0023 p-v*lua= 0.057 



pare ant rad contrast 

ncn-otoOqu* a* 0.1722 b= 0.02344 obtrqua a* 020257 b= 0.03495 
thats- -1 7.602 x* -0.0054 y- 0.0093 p-valua- 0268 



Figure 4.4. Ellipses and values for Subjects 1-4 at 3 cpd 



Ellipse Type 


Subj 


a 


b 


theta 


centers 


center.y 


delta.a 


delta.b 


*p-value 


non-oblique 


1 


0.05594 


0.02284 


-0.3235 


0.01224 


0.0001 








oblique 


1 


0.06445 


0.02256 


-0.3235 


0.01224 


0.0001 


0.0085 


-0.0003 


0.028 


non-oblique 


2 


0.0459 


0.01365 


-0.2959 


-0.0153 


0.01315 








oblique 


2 


0.0556 


0.01554 


-0.2959 


-0.0153 


0.01315 


0.0097 


0.0019 


0.005 


non-oblique 


3 


0.08941 


0.02025 


-0.2149 


0.00796 


0.00232 








oblique 


i 3 


0.10976 


0.02385 


-0.2149 


0.00796 


0.00232 


0.0203 


0.0036 


0.057 


non-oblique 


4 


0.1722 


0.02344 


-0.3072 


-0.0054 


0.0093 








oblique 


4 


0.14182 


0.03495 


-0.3072 


-0.0054 


0.0093 


-0.0304 


-0.0584 


0.268 


*Ho delta. a=0 




















Ha delta.a<>0 


i 

















Table 4.3. Ellipses and values for Subjects 1-4 at 3 cpd 
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p*rc*ntgr#*n contrast 



Sub 1 7 cpd 2 h z 

non-oblique = points/solid line oblique = pluses/dashed line 




norvcWiqua m- 0.06371 b- 0.02381 oWfqua •- 0.06646 b= 0.03584 
lhata- -20.009 x» 0.0052 y> 0.0001 p-vaJua- 0.588 



Sub 2 7 cpd 2 h z 



non-oblique = points/solid line oblique = pluses/dashed line 




parcant rad contrast 

non-ctttqua 007852 b* 0,01751 obttqua a> 0.09241 b“ 0.02723 
that!" -16 684 x« -0.006 y* 0.008 p-valua* 0.087 



Sub 3 7 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 




Sub 4 7 cpd 2 hz 

non-oblique = points/solid line oblique = pluses/dashed line 




pare ant rad contrast 

non-oWiqua a* 5 .35208 b* 0.0178 obliqua a« 0.11977 b» 0.02158 
thata- -14.176 x» 4) 0077 y» 0.0024 p-vaJua» 0.059 



parcant rad contrast 

non-obiqua a« 021246 b= 0.04223 obfiqua a- 0.1788 b* 0.04969 
thata* -1 3.456 x» 0.0081 y- 0.0018 p-valua- 0.43 



Figure 4.5. Ellipses and values for Subjects 1-4 at 7 cpd 



Ellipse Type 


Subj 


a 


b 


theta 


center jc 


center.y 


delta.a 


delta.b 


*p-value 


non-oblique 


1 


0.06371 


0.02391 


-0.3492 


0.00517 


0.00012 








oblique 


1 


0.06646 


0.03584 


-0.3492 


0.00517 


0.00012 


0.0028 


0.0119 


0.588 


non-oblique 


2 


0.07852 


0.01751 


-0.2912 


-0.006 


0.00799 








oblique 


2 


0.09241 


0.02723 


-0.2912 


-0.006 


0.00799 


0.0139 


-0.0447 


0.087 


non-oblique 


3 


5.35208 


0.0178 


-0.2474 


-0.0077 


0.00238 








oblique 


3 


0.11978 


0.02158 


-0.2474 


-0.0077 


0.00238 


-5.2323 


-0.0394 


0.059 


non -oblique 


4 


0.21246 


0.04223 


-0.2349 


0.00806 


0.00179 








oblique 


4 


0.1788 


0.04969 


-0.2349 


0.00806 


0.00179 


-0.0337 


0.0075 


0.43 


*Ho delta. a=0 




















Ha delta. a<>0 





















Table 4.4. Ellipses and values for Subjects 1-4 at 7 cpd 
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1 cpd 


3 cpd 


7 cpd 




run 


non-oblique 


oblique 


non-oblique 


oblique 


non-oblique 


oblique 


Subject 1 


i 


0.359 


0.436 


0.000 


0.562 


0.000 


0.397 




2 


0.056 


0.199 


0.242 


0.048 


0.191 


0.393 




3 


0.818 


0.544 


0.680 


0.000 


0.558 


0.933 




4 


0.399 


0.673 


0.152 


0.181 


0.529 


0.001 




5 


0.066 


0.831 


0.722 


0.141 


0.762 


0.975 


















Subject 2 


1 


0.065 


0.695 


0.812 


0.758 


0.000 


0.000 




2 


0.070 


0.800 


0.005 


0.620 


0.899 


0.146 




3 


0.222 


0.019 


0.900 


0.014 


0.462 


0.013 




4 


0.002 


0.004 


0.002 


0.676 


0.000 


0.126 


















Subject 3 


1 


0.000 


0.604 


0.021 


0.759 


0.834 


0.217 




2 


0.770 


0.033 


0.796 


0.095 


0.006 


0.231 




3 


0.226 


0.017 


0.079 


0.001 


0.101 


0.004 


















Subject 4 


1 


0.320 


0.752 


0.057 


0.394 


0.000 


0.071 




2 


0.777 


0.258 


0.546 


0.042 


0.513 


0.000 




3 


0.397 


0.926 


0.415 


0.000 


0.136 


0.000 




4 


0.620 


n/a 


n/a 


0.659 


0.003 


n/a 



Table 4.5. P-values for comparing one subjects run against their other runs at that same condition. 
P-values have not been corrected for multiple comparisons. Bold values are less than ,05/(number of 
runs at that condition). 
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V. CONCLUSIONS 



The purpose of this thesis was to investigate the human visual system’s ability to 
detect certain simple targets. This thesis investigated how an object’s spatial, temporal, and 
color features affected humans’ detection of objects. The results showed that certain spatial 
and chromatic qualities do, indeed, inhibit detection. While real-world objects are much more 
complex than laboratory stimuli, knowledge of spatial and chromatic qualities that inhibit 
detection will assist military designers in the quest for better camouflage. 

Other studies of use to military designers include the numerous studies documenting 
an achromatic Class-I oblique effect and the fact that it is generally found psychophysically 
only at high spatial frequencies. This study produced similar results with all subjects 
displaying an achromatic Class-I oblique effect (p-values of 0, 0, .139 and 0) at a spatial 
frequency of 7 cpd. Previous studies documenting a chromatic Class-I oblique effect, or lack 
thereof, are less useful due to conflicting results and possible problems with luminance 
artifacts tainting results (Kelly, 1976; Murasagi and Cavanagh, 1988). Indeed, the work done 
by Kelly and Murasagi and Cavanagh highlighted the problems in determining a chromatic 
oblique effect due to the difficulty of obtaining isoluminance for a subject. Any deviation can 
lead to the introduction of luminance artifacts and can corrupt the results of the experiment. 
The methodology used in this thesis takes advantage of the elliptical shape of the curve 
connecting thresholds at a fixed temporal and spatial frequency, and makes the exact 
determination of isoluminance unnecessary. 

This thesis supports the hypothesis that a Class-I chromatic oblique does exist. At a 
spatial frequency of 3 cpd, a chromatic oblique effect is evident. A chromatic oblique effect 
is shown for three of the four subjects and, with an alpha of .05, is significant for two of the 
four subjects. Additionally, while the p-value for Subject Three (.057) is not less than .05, 
this subject conducted only three runs with 12 thresholds. An additional run would likely 
reduce this p-value to a value less than .05. 
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The main value of this study is the tool it provides for further investigation of a Class-I 
chromatic oblique effect without the problems associated with luminance artifacts. However, 
this tool does have its drawbacks. The data collection required is extremely time-intensive. 
A total of at least five runs at each spatial frequency is desirable due to the variability of the 
data. With each session lasting approximately 45 minutes to an hour and one run consisting 
of an oblique session and a non-oblique session, the time required for five runs is 
approximately 7.5 to ten hours. This time does not account for three sessions needed to 
determine a subject’s maximum oblique and non-oblique sensitivity. This large time 
commitment on the part of subjects poses problems, as their motivation begins to wane. 
Motivation was a possible problem, as UL subjects did not run consistently; their average run 
time (days) was 19 and 18 versus 10 and 12 for the NPS subjects. This undoubtedly affected 
the variability of their data, as evidenced by the number of their p-values in Table 4.5 less than 
.05. 

While the determination of exact isoluminance is not required, it is desirable to 
determine an approximate isoluminance axis. If an ellipse major to minor axis ratio is four 
or greater, then Sellers et al. (1986) state that the major axis is a good approximation of the 
equiluminous axis. Major to minor axis ratios in this study averaged 3.3 and exceeded four 
only one-third of the time. For Subject One at 7 cpd graphically (Figure 4.3), it appears that 
a Class-I chromatic oblique effect is evident. Taking into account that the highest major to 
minor axis ratio for this subject at this spatial frequency is 2.66, a chromatic oblique effect is 
even more likely since the equiluminous axis is not likely approximated that well by the major 
axis. In this case, the p-value of .588 does not provide much information regarding what 
actually occurred. The statistical test resulting from the model formulation tests the 
significance in the difference of the length of the axes and does not account for the fact that 
the true isoluminance axis may not be aligned with the major axis. 

This thesis has provided an excellent tool for further research. Possible improvements 
include determination of a subject’s equiluminous axis prior to running the experiment, thus 
enabling the experimenter to choose ratios that would run through this approximate axis. A 
spatial frequency of 3 cpd is worthy of further study with more runs and fewer confounding 
variables (different monitors). 
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APPENDIX A. S-PLUS CODE 



> ellipse 

function(x, y, a = 2, b = 3, e = 0, theta = 0, fit. center = F, grad = T, 
is . there . hess = T, plot.it = T, chat = T) 

{ 

# This is for a Class-I type ellipse 

# Fit a least-squares ellipse centered at 0 with semi-axes (a, b) 

# and angle to the origin theta, to the data in x, y. The ellipse 

# is here parameterized by a, e (the eccentricity) and theta, 

# in that order, a is always > b. 

# a,b,e and theta are starting points 

# fit.center=F pins the ellipse to the origin when this is true the 

# ellipse is allowed to "float" 

# grad=gradient 

# . hess=hessian 

# plot.it activates plot 

# chat=T shows values as they are computed 

# 

# If a is supplied, and it 1 s a vector, then we've been given 

# starting points for all the parameters. Use 'em, first making 

# sure that there is the right number (3 if we're not fitting 

# the center, and 5 if we are) . 

# 



if ( !missing (a) && length(a) > 1) { 

if (( fit . center && length(a) != 5) || (! fit. center && 

length (a) ! = 3 ) ) 

stop (paste ( "Parameter vector has length", length (a), 
", expecting ", if else ( fit . center, 5, 3))) 
if (length (names (a [2] ) ) — 0 | | names (a [2]) == "e") { 

if (length (names (a [2] ) ) == 0) 

warning ( "No param names: using e in pos. 2") 
e <- a [ 2 ] 

b <- a[l] * sqrt(l - e A 2) 

} 

else { 

b <- a [2] 

e <- sqrt (1 - (b/a[l] ) A 2) 

} 

theta <- a[3] 
if ( fit . center) { 

center. x <- a[4] 
center. y <- a[5] 

} 

a <- a[l] 



else { 

if (missing (a) ) 

a <- 0.5 * diff (range (x) ) 
if (missing (b) ) 

b <- 0.5 * diff (range (y) ) 

e <- ifelse(a > b, sqrt(l - (b/a) A 2), sqrt(l - (a/b) A 2)) 
if (missing (theta) ) { 

Is. out <- lsfit(x, y) 

theta <- atan(ls.out$coef [2] ) 
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} 

center. x <- mean(x) 
center. y <- mean(y) 

} 

tt <- ell . tt (x, y) 
if(plot.it) { 

graph <- dev. list ()[ "win. graph" ][ 1] 
if (length (graph) == 0 I | all (is . na (graph )) ) 
win . graph ( ) 
else dev. set (graph) 



if (grad) 

grad.func <- ell. grad 
else grad.func <- NULL 
if ( fit . center) { 

start. vec <- c (a = a, e = e, theta = theta, center. x = 
center. x, center. y = center. y) 
lower. vec <- c(0.0001, 0.0001, -2 * pi, - Inf, - Inf) 
upper. vec <- c (Inf , 0.999999, 2 * pi, Inf, Inf) 



else { 

start. vec <- c (a = a, e = e, theta = theta) 
lower. vec <- c(0.0001, 0.0001, -2 * pi) 

upper. vec <- c(Inf, 0.999999, 2 * pi) 

} 

out <- nlminb (start = start. vec, objective = ell. res, gradient = 

grad.func, hessian = is . there . hess , lower = lower. vec, upper 
= upper. vec, tt = tt, my.x = x, my. y = y, plot.it = 
plot.it, chat = chat, is. there. h ess = is . there . hess, 
fit. center = fit. center, step. min = 100 * 
.Machine3double.eps, scale. upd =1) # 

### p. names <- names (out $parameters) 

### cat ("In ellipse, check p.names\n") 

### browser () 

b <- out$parameters [ "a"] * sqrt(l - out$parameters[ ,, e"] A 2) 
if (length (out$parameters) > 3) { 

out$parameters <- c (out$parameters [ "a" ] , b = b, out$ 
parameters [ 3 : length (out$parameters ) ] ) # 

names (out$parameters) <- c("a", "b", "theta", "center. x", 

"center. y") # Beats me... 

} 

else { 

out$parameters <- c (out$parameters [ "a"] , b = b, out$ 
parameters [ 3] ) # 

names (out$parameters) <- c("a", "b", "theta") 

} 

out$sigma.sq <- out$obj /length (x) 
out$sigma <- sqrt (out$sigma . sq) 
out$aux <- NULL 
out$x <- x 
out$y <- y 

if (length ( out$hessian) > 0) { 

if (qr (out$hessian) $rank < ncol (out$hessian) ) 
out$cov <- "Can't invert Hessian" 
else out$cov <- out$sigma.sq * solve (out$hessian) 



out$tt <- tt 

out$fitted. tt <- ell.tt(x - center. x, y - center. y) 
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=#= 4*= :#= 41= 4*= =H= 4t= = 8 = 4t= 4*= ^ 



pred <- ell.pred(out$fitted. tt, out$parameters [ "a"] , 
out$parameters[ ,, b"] , out$parameters [ " theta" ] , 
return. unrotated. too = F, fit . center = fit. center, 
center. x = if else ( fit . center, out$ 

parameters [ "center. x"] , 0), center. y = ifelse ( fit . center , 
out$ parameters [ "center. y " ] , 0)) 

out$fitted.x <- pred$x 
out$fitted.y <- pred$y 

out$fitted.r <- sqrt (pred$x A 2 + pred$y A 2) 
class (out) <- "ellipse" 
return (out) 

} 



> ellipse . II 

function (x, y, a = 2, b = 3, theta = 0, delta = 0, fit .-center = F, grad 
= T, is. there. hess = T, plot.it = T, chat = T, class.I = rep(T, 
length (x) ) , which. type = 1) 

{ 

This is for a Class-II type ellipse 

class.I seperates oblique and non-oblique data. For an x vector of 
length 160 where the first 80 points were non-oblique data the 
class.I vector should consist of a boolean vector or length 160 # 
comprised of 80 T's followed by 80 F' s 
which. type=l when testing differences in the major axes (a's or # 
chromaticity) 

which. type=2 when testing differences in the minor axes (b' s 
or luminance) 

This version of ellipse works, but you must set grad=F, is . there . hess=F 
and plot.it=F 

Fit a least-squares ellipse centered at 0 with semi-axes (a, b) 
and angle to the origin theta, to the data in x, y. The ellipse 
is now parameterized by a, b (not e) and theta, in that order, 
a is always > b; we can enforce that at the end. 

if (is .matrix (x) && ncol(x) > 1) { 

if (any (dimnames (x) [ [2] ] == "y") ) { 

y <- x [ / "y"] 

if (any (dimnames (x) [ [2] ] == "x") ) 
x <— X [ , "x"] 

else x <- x[, 1] 

} 

else { 

y <- x[, 2] 
x <- x[, 1] 

} 

) 

if (is . list (x) ) { 

if (any (names (x) == "y") ) { 
y <- x$y 

if (any (names (x) == "x")) 
x <- x$x 
else x <- x[ [1] ] 

} 

else { 

y <- x[[2]] 

X <- x[[l]] 
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} 



} 



# 

# 

# 

# 

# 

# 



If a is supplied, and it's a vector, then we've been given 
starting points for all the parameters. Use 'em, first making 
sure that there is the right number (3 if we're not fitting 
the center, and 5 if we are) . These #'s increase 1 for every delta 
estimated. 

if ( Imissing (a) && length(a) > 1) { 

if (( fit . center && length (a) != 5) || (! fit-center && 

length (a) != 3) ) 

stop (paste ( "Parameter vector has length", length (a), 
", expecting ", i f else ( fit . center , 5, 3))) 

b <- a [2 ] 

e <- sqrt (1 - (b/a[l] ) ^2) 

theta <- a [3] 
if ( fit . center) { 

center. x <- a[4] 
center. y <- a[5] 

} 

a <- a[l] 



} 

else { 

if (missing (a ) ) 

a <- 0.5 * dif f ( range (x) ) 
if (missing (b) ) 

b <- 0.5 * dif f (range (y) ) 
if ( a < b) { 

switcheroo o a 
a <- b 

b <- switcheroo 



} 

e <- sqrt(l - (b/a)~2) 
if (missing ( theta ) ) { 

Is. out <- lsfit(x, y) 

theta <- atan (Is . out$coef [2 ] ) 



center. x <- mean(x) 
center. y <- mean(y) 

} 

tt <- ell.tt(x, y) 
if (plot . it) { 

graph <- dev.list() [ "win . graph"] [1] 
if (length (graph) == 0 | | all ( is . na (graph) ) ) 
win . graph ( ) 
else dev. set (graph) 

} 

if (grad) 

grad.func <- ell. grad. II 
else grad.func <- NULL 
if ( fit . center) { 

start. vec <- c(a = a, b = b, theta = theta, center. x = 
center. x, center. y = center. y, delta = delta) 
lower. vec <- c(0. 00001, 0.00001, -2 * pi, - Inf, - Inf, 
Inf) 

upper. vec <- c(Inf, Inf, 2 * pi. Inf, Inf, Inf) 

} 

else { 
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start. vec <- c(a = a, b = b, theta = theta, delta = delta) 
lower. vec <- c(0. 00001, 0.00001, -2 * pi, - Inf) 

upper. vec <- c(Inf, Inf, 2 * pi. Inf) 

} 

out <- nlminb (start = start. vec, objective = ell. res. II, gradient 
= grad.func, hessian = is . there .hess, lower = 
lower. vec, upper = upper. vec, tt = tt, my.x = x, my. y 
= y, plot.it = plot.it, chat = chat, is . there . hess = 
is. there. h ess, fit. center = fit. center, class.I = 
class. I, which. type = which. type, step. min = 100 * 
,Machine$double . eps, scale. upd = 1) # 

### p. names <- names (out$parameters ) 

### cat ("In ellipse, check p.names\n") 

### browser () 

if (length (out$parameters) > 4) 

names ( out$parameters ) <- c("a", "b", "theta", "center. x", 

"center. y", "delta") # Beats me... 
else names (out$parameters) <- c("a", "b", "theta", "delta") 

out$sigma.sq <- out$obj /length (x) 
out$sigma <- sqrt ( out$sigma . sq) 
out$aux <- NULL 
out$x <- x 
out$y <- y 

if (length (out$hessian) > 0) { 

if (qr (out$hessian) $rank < ncol ( out$hessian) ) 
out$cov <- "Can’t invert Hessian" 
else out$cov <- out$sigma . sq * solve ( out$hessian) 

} 

out$tt <- tt 

out$fitted. tt <- ell.tt(x - center. x, y - center. y) 
pred<-ell .pred (out$fitted. tt, out$parameters ["a"] , 
out$parameters [ "b"] , out$parameters [ "theta"] , 
return. unrotated. too = F, fit. center = fit. center, 
center. x = if else ( fit . center , out$parameters [ "center . x" ] , 

0), center. y = if else ( fit . center, out$parameters [ "center . y" ] , 
0)) 

out$fitted.x <- pred$x 
out$fitted.y <- pred$y 

out$fitted.r <- sqrt (pred$x^2 + pred$y / '2) 
class (out) <- "ellipse" 
return (out ) 



> ellipse. Ill 

function(x, y, a = 2, b = 3, theta = 0, delta. a = 0, delta. b = 0, 

fit. center = F, grad = T, is . there . hess = T, plot.it = T, chat = 

T, class.I = rep(T, length(x))) 

{ 

# This is for Class-Ill ellipses 

# This version of ellipse works, but you must set grad=F, is . there . hess=F 

# and plot.it=F 

# 

# 

# Fit a least-squares ellipse centered at 0 with semi-axes (a, b) 

# and angle to the origin theta, to the data in x, y. The ellipse 

# is now parameterized by a, b (not e) and theta, in that order. 

# a is always > b; we can enforce that at the end. 
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# 



# 

# 

# 

# 

# 

# 



if (is .matrix (x) && ncol (x) > 1) { 

if (any (dimnames (x) [ [2] ] == "y")) { 

y <- x[, "y"] 

if (any (dimnames (x) [ [2] ] == "x")) 
x <- x[, n x n ] 
else x <- x[, 1] 

} 

else { 

y <- x[, 2] 
x <- x[, 1] 

} 

} 

if (is . list (x) ) { 

if (any (names (x) == "y") ) { 

y <- x$y 

if (any (names (x) == "x" ) ) 
x <- x$x 
else x <- x [ [ 1] ] 

} 

else { 

y <- x[[2]] 

X <- x[[l] ] 

} 

} 



If a is supplied, and it's a vector, then we've been given 
starting points for all the parameters. Use * em, first making 
sure that there is the right number (3 if we * re not fitting 
the center, and 5 if we are) . These #'s increase 1 for every delta 
estimated 

if ( Imissing (a) && length(a) > 1) { 

if (( fit . center && length (a) != 5) || (! fit. center && 

length (a ) != 3) ) 

stop (paste ( "Parameter vector has length", length (a), 
", expecting ", if else ( fit . center , 5, 3))) 

b <- a [2 ] 

e <- sqrt ( 1 - (b/a[l]) A 2) 

theta <- a [3] 
if ( fit . center ) { 

center. x <- a[4] 
center. y <- a[5] 

} 

a <- a [ 1 ] 



} 

else { 

if (missing (a) ) 

a <- 0.5 * diff (range (x) ) 
if (missing (b ) ) 

b <- 0.5 * diff (range (y) ) 
if (a < b) { 

switcheroo <- a 
a <- b 

b <- switcheroo 



} 

e <- sqrt (1 - (b/a) A 2) 

if (missing (theta) ) { 

Is . out <- lsfit(x, y) 



50 



theta <- atan (Is . out$coef [2] ) 



} 

center. x <- mean(x) 
center. y <- mean(y) 

) 

tt <- ell.tt(x, y) 
if(plot.it) { 

graph <- dev. list ()[ "win. graph"] [ 1] 
if (length (graph) == 0 | | all (is . na (graph) ) ) 
win. graph ( ) 
else dev. set (graph) 

) 

if (grad) 

grad.func <- ell. grad. Ill 
else grad.func <- NULL 
if ( fit . center) { 

start. vec <- c(a = a, b = b, theta = theta, center. x = 
center. x, center. y = center. y, delta. a = delta. a, 
delta. b = delta. b) 

lower. vec <- c(0. 00001, 0.00001, -2 * pi, - Inf, - Inf, 0, 
0) 

upper. vec <- c(Inf, Inf, 2 * pi. Inf, Inf, Inf, Inf) 

) 

else { 

start. vec <- c ( a = a, b = b, theta = theta, delta. a = 
delta . a, delta . b = delta. b) 
lower. vec <- c(0. 00001, 0.00001, -2 * pi, 0, 0) 

upper. vec <- c(Inf, Inf, 2 * pi. Inf, Inf) 

) 

out <- nlminb (start = start. vec, objective = ell. res. Ill, gradient 
= grad.func, hessian = is . there . hess, lower = lower. vec, 
upper = upper. vec, tt = tt, my.x = x, my.y = y, plot.it = 
plot.it, chat = chat, is . there . hess = is . there . hess, 
fit. center = fit. center, class. I = class. I, step. min = 100 * 
. Machine$double . eps, scale. upd = 1) # 

### p. names <- names (out$parameters) 

### cat ("In ellipse, check p.names\n") 

### browser () 

if (length (out$parameters) > 4) 

names ( out$parameters) <- c("a", "b", "theta", "center. x", 

"center. y", "delta. a", "delta. b") # Beats me... 

else names (out$parameters ) <- c("a", "b", "theta", "delta. a", 

"delta. b") 

out$sigma.sq <- out$obj /length (x) 
out$sigma <- sqrt (out$sigma . sq) 
out$aux <- NULL 
out$x <- x 
out$y <- y 

if ( length (out$hessian) > 0) { 

if ( qr (out$hessian) $rank < ncol (out$hessian) ) 
out$cov <- "Can't invert Hessian" 
else out$cov <- out$sigma.sq * solve ( out$hessian) 

} 

out$tt <- tt 

out$fitted. tt <- ell.tt(x - center. x, y - center. y) # 
pred <- ell .pred (out$fitted. tt, out$parameters [ "a"] , 
out$parameters [ 

"b" ], out$parameters [ "theta" ] , return . unrotated . too = F, 



### 

### 

### 
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### fit. center = fit. center, center. x = if else ( fit . center , out$ 

### parameters [ "center . x"] , 0), center. y = if else ( fit . center , 

### out$parameters [ "center . y"] , 0) ) 

### out$fitted.x <- pred$x 

### out$fitted.y <- pred$y 

### out$fitted.r <- sqrt (pred$x^2 + pred$y^2) 

class (out) <- "ellipse" 
return (out ) 

} 



> ell . res 

function (params, tt, my.x, my.y, is . there . hess, fit. center, plot.it, 
chat) 

{ 

# 

# ell. res: Compute objective to be minimized. 

# 

# This computes the objective function: the sum of squared 

# differences between the observed points on the ellipse 

# (after transformation) and the predicted ones. 

# 

# "params" is the vector (a, e, theta) . Get them out, and 

# compute rat, the ratio a/b. 

# 

a <- params [1] 
e <- params [2] 
if (e > 0.99) 

return (1000) 
b <- a * sqrt ( 1 - e A 2) 
theta <- params [3] # 

if ( fit . center == T) { 

center. x <- params [4] 
center. y <— params [5] # 

tt <- ell. tt (my.x - center. x, my.y - center. y) 



else { 

center. x <- center. y <- 0 

) 

fitted. r <- ell.pred(tt, a, b, theta, fit. center = fit. center, 
center. x = center. x, center. y = center. y) # 
new.x <- fitted. r$x 
new.y <- fitted. r$y # 

# Plot it: add dotted lines at x = 0 and y = 0. 

# 

if(plot.it) { 

plot (my.x, my.y, xlim = range (my.x, new.x), ylim = 
range(my.y, new.y)) 
abline(h = 0, lty = 2) 
abline ( v = 0, lty =2) # 

points (new. x, new.y, pch = 1, col =4) # 

points (center . x, center. y, pch = 1, col = 2) 

} 

### ford in 1 : length (tt) ) 

### polar (c ( 0, 0.5), rep(tt[i], 2), type = "1") # 

# 

# Get fitted x and y; compute and return objective. 

# 
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obj <- sum( (my.x - new.x) A 2) + sum((my.y - new.y) A 2) # 
if (chat) cat ( "a : ", signif (a, 4), "/b:", signif (b, 4), ",th:", 

signif (theta, 4), if else ( fit . center , paste ("; x, y: ", signif ( 
center. x, 4), signif (center . y, 4)), "") , "/obj:", 

signif (obj, 4), "\n") ### cat("BTW: ") 

### ell . grad (params, tt, my.x, my.y, is . there . hess ) 

### cat ( "\n" ) 

return (obj ) 

} 

> ell . res . II 

function (params, tt, my.x, my.y, is . there . hess, fit. center, plot.it, 
chat, class. I, which. type) 

{ 



ell. res: Compute objective to be minimized. This version is the 
Class-II one. 

This computes the objective function: the sum of squared 
differences between the observed points on the ellipse 
(after transformation) and the predicted ones. 

"params" is the vector (a, b, theta) . 

a <- params [1] 

b <- params [2] 

theta <- params [3] # 

if (is .null (class . I ) ) { 

class. I <- rep(T, length (my . x) ) 
delta <- 0 

} 

if (fit. center == T) { 

center. x <- params [4] 
center. y <- params [5] # 

delta <- params [6] 

tt <- ell. tt (my.x - center. x, my.y - center. y) 

} 

else { 

center. x <- center. y <- 0 
delta <- params [4] 

} 

if ( sum (class . I ) < length (my . x) ) { 

fitted, r. I <*- ell .pred (tt [class . I] , a, b, theta, fit. center 
= fit. center, center. x = center. x, center. y = 
center. y) # 
if (which . type == 1) 

fitted. r. II <- ell .pred (tt [! class . I] , a + delta, b, 
theta, fit. center = fit. center, center. x = 
center. x, center. y = center. y) 
else fitted. r. II <- ell .pred ( tt [ I class . I] , a, b + delta, 
theta, fit. center = fit. center, center. x = 
center. x, center. y = center. y) # 

} 

else fitted. r. I <- ell.pred(tt, a, b, theta, fit. center = 
fit. center, center. x = center. x, center. y = center. y) # 
new.x <- numeric (length (my . y) ) 
new.y <- numeric (length (my . y) ) 
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new. x [class . I ] <- fitted. r.I$x 
new. y [class . I] <- fitted. r.I$y 
if (sum( ! class . I) > 0) { 

new. x [ ! class . I ] <- f itted . r . II$x 
new. y [ ! class . I] <- f itted . r . II$y 



If plot it, add dotted lines at x = 0 and y = 0, plus points. 
if(plot.it) { 

plot (my. x, my.y, xlim = range (my. x, new.x), ylim = 
range (my. y, new.y) ) 

abline(h = 0, lty = 2J 
abline(v = 0, lty =2) # 

points (new. x, new.y, pch = 1, col =4) # 

points (center . x, center. y, pch = 1, col = 2) 

} 



cat ( "grad. norm is ", sum (ell . grad. II (params, tt, my.x, my.y, 

is . there . hess , fit. center, class. I, which . type) A 2 ) , "\n") 



Get fitted x and y; compute and return objective. 

obj <- sum ( (my.x - new.x) A 2) + sum((my.y - new.y) A 2) # 

if (chat) 

cat ("a:", signif(a, 4), ", delta: ", signif (delta, 4), 

", b: ", signif (b, 4), ",th:", signif (theta, 4), ifelse( 

fit. center, paste ("; x, y :" , signif (center . x, 4), 
signif (center . y, 4)), ""), ";obj:", signif (obj, 4), 
"\n" ) 

return (obj ) 



ell. res. Ill 
unction (params, tt, 
chat, class.I) 



my.x, my.y. 



is .there. hess. 



fit. center, plot.it. 



ell. res: Compute objective to be minimized. This version is the 
Class-Ill one. 

This computes the objective function: the sum of squared 
differences between the observed points on the ellipse 
(after transformation) and the predicted ones. 

"params" is the vector (a, b, theta) . 

a <- params [1] 

b <- params [2] 

theta <- params [3] # 

if (is .null (class . I ) ) { 

class.I <- rep(T, length (my . x ) ) 
delta. a <- delta. b <- 0 

} 

if ( fit . center == T) { 
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center. x <- params[4] 
center. y <- params[5] # 
delta. a <- params[6] 
delta.b <- params[7] 

tt <- ell.tt(my.x - center.x, my.y - center. y) 

} 

else { 

center.x <- center. y <- 0 
delta. a <- params[4] 
delta.b <- params[5] 

) 

if (sum (class . I ) < length (my . x) } { 

fitted. r. I <- ell . pred ( tt [class . I] , a, b, theta, fit. center 
= fit. center, center.x = center.x, center. y = 
center. y) # 

fitted. r. II <- ell . pred ( tt [! class . I] , a + delta. a, b + 

delta.b, theta, fit. center = fit. center, center.x = 
center.x, center. y = center. y) # 

} 

else fitted. r. I <- ell.predftt, a, b, theta, fit. center = 

fit. center, center.x = center.x, center. y = center. y) # 
new.x <- numeric (length (my . y) ) 
new.y <- numeric (length (my . y) ) 
new.x [class . I] <- fitted. r.I$x 
new. y [class . I ] <- fitted. r.I$y 
if ( sum ( ! class . I ) > 0) { 

new. x [ ! class . I ] <- fitted. r . II$x 
new. y [ ! class . I] <- fitted. r . II$y 

) 

# 

# 

# If plot it, add dotted lines at x = 0 and y = 0, plus points. 

# 

if (plot . it) { 

plot (my. x, my.y, xlim = range (my. x, new.x), ylim = 
range (my.y, new.y) ) 
abline(h = 0, lty = 2) 
abline (v = 0, lty =2) # 

points (new. x, new.y, pch = 1, col =4) # 

points (center . x, center. y, pch = 1, col = 2) 

} 

### 

## cat ("grad. norm is ", sum ( ell . grad. II (params, tt, my.x, my.y, 

## is . there . hess, fit. center, class.I, which . type) A 2 ) , "\n") 

# 

# Get fitted x and y; compute and return objective. 

# 

obj <- sum ( (my.x - new.x) A 2) + sum((my.y - new.y) A 2) # 
if ( chat) 

cat ("a:", signif(a, 4), ", delta. a: ", signif (delta . a, 4), 

",b:", signif (b, 4), "delta.b: ", signif (delta . b , 4), 
",th:", signif (theta, 4), if else ( fit . center , paste ( 
";x,y:", signif (center . x, 4), signif (center . y, 4)), 
""), "; obj : " , signif (obj, 4), "\n") 

return (obj ) 

} 
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> ell . grad 

function (params, tt, my.x, my.y, is . there . hess , fit. center) 

{ 

tol <- le-006 
a <- params [1] 
e <- params [2] 
theta <- params [3] 
if (fit . center == T) { 

center. x <- params [4] 
center. y <- params [5] 

tt <- ell. tt (my.x - center. x, my.y - center. y) 

} 

else { 

center. x <- center. y <- 0 



### 

### 

### 



(tt 

(tt 



theta) ) 
theta) ) 



'2 

'2 



b <- a * sqrt(l - e A 2) 

fitted <- ell.pred(tt, a, b, theta, return . unrotated. too = T, 
fit. center = fit. center, center. x = center. x, center. y = 
center . y) 

xprime <- fitted$x . prime 
yprime <- fitted$y . prime 
x <- fitted$x 
y <- fitted$y 
cos. theta <- cos (theta) 
sin. theta <- sin (theta) 
cos . 2 . tt . theta <- cos (2 
sin . 2 . tt . theta <- sin (2 
sinsq. tt. theta <- (sin(tt - theta) ) A 2 
cossq. tt . theta <- (cos (tt - theta) ) A 2 
sinsq. 2 . tt . theta <- (sin (2 * (tt - theta))) 
consq. 2 . tt . theta <- (cos (2 * (tt - theta))) 
one . minus . e . sq <- 1 - e A 2 

denom <- cossq. tt . theta * one . minus . e . sq + sinsq . tt. theta 
dxprime.da <- xprime/a 

dxprime.de <- ( ( - a A 2/(4 * xprime)) * (e * 
sins q. 2. tt. theta) ) / (denom A 2) 
dxprime . de [abs (xprime) < tol] <- 0 
dyprime.da <- yprime/a 
dyprime.de <- - (one .minus . e . sq 

-xprime A 2) ) /yprime 
dyprime . de [abs (yprime) < tol] <- 
dx.da <- cos. theta * dxprime.da 

dx. de <- cos. theta * dxprime.de 

dy. da <- sin. theta 
dy.de <- sin. theta 

x. diff <- my.x - x 

y. diff <- my.y - y # 

grad. mat <- matrix (0, length (x) 
grad.mat[, 1] <- -2 * (x.diff * 
grad.mat[, 2] <- -2 * (x.diff * 
grad. a <- -2 * sum(x.diff * dx.da + y 
grad.e <- -2 * sum(x.diff * dx.de + y 
num <- one . minus . e . sq * sin(2 * (tt - 
dxprime . dtheta <- (a A 2/ (2 * xprime)) 

dxprime . dtheta [abs (xprime ) < tol] <- 
dyprime . dtheta <- - ( one . minus . e . sq 

yprime 

dyprime . dtheta [abs (yprime) < tol] <- 



xprime * dxprime.de + e * (a A 2 



* dxprime.da 

* dxprime.de 



sin . theta 
sin . theta 
cos . theta 
cos . theta 



dyprime . da 
dyprime . de 
dyprime . da 
dyprime . de 



3) 

dx.da + y.diff * dy.da) 
dx.de + y.diff * dy.de) 
diff * dy.da) 
diff * dy.de) 
theta) ) 
(num/denom A 2 ) 



xprime * dxprime . dtheta ) / 
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### 

### 

### 



dx. dtheta <- - (y - center. y) + cos. theta * dxprime . dtheta - 

sin. theta * dyprime . dtheta 

dy. dtheta <- (x - center. x) + sin. theta * dxprime . dtheta + 
cos. theta * dyprime . dtheta 

grad. theta <- -2 * sum(x.diff * dx.dtheta + y.diff * dy. dtheta) 

# 



grad.mat[, 3] <- -2 * (x.diff * dx.dtheta + y.diff * dy. dtheta) 
cat ("Grad mat approx. is...\n" ) 
print (t (grad. mat ) %*% grad. mat) 
if (fit . center == F) 

grad <- c (grad. a, grad.e, grad. theta) 

else { 



dxprime. dt <- - dxprime . dtheta 

dyprime. dt <- - dyprime . dtheta 

R.sq <- (my.x - center. x) A 2 + (my.y - 
dt.dxO <- (my.y - center . y) /R. sq 
dt.dyO <- - (my.x - center . x) /R. sq 

dxprime. dxO <- dxprime. dt * dt.dxO 
dxprime. dyO <- dxprime. dt * dt.dyO 
dyprime. dxO <- dyprime.dt * dt.dxO 
dyprime. dyO <- dyprime.dt * dt.dyO # 
dyprime. dxprime <- one .minus . e . sq * ( 
dyprime. dxprime [abs (yp rime) < tol] <- 

dx. dxO <- (cos. theta * dxprime. dxO) - 

dyprime.dxO) + 1 

dy. dxO <- (sin. theta * dxprime. dxO) + 

dyprime . dxO) 

dx. dyO <- (cos. theta * dxprime. dyO) - 

dyprime. dyO) 

dy. dyO <- (sin. theta * dxprime. dyO) + 

dyprime. dyO) + 1 

grad.xO <- -2 * sum (x.diff * dx.dxO + 
grad.yO <- -2 * sum(x.diff * dx.dyO + 
grad <- c (grad. a, grad.e, grad. theta. 



center . y) A 2 



- xp rime/ yp rime) 

0 # 

(sin. theta * 

(cos. theta * 

(sin. theta * 

(cos. theta * 

y.diff * dy.dxO) 
y.diff * dy.dyO) 
grad.xO, grad.yO) 



} 

if (is . there . hess == F) 
return (grad) 

d2xprime.da2 <- d2yprime.da2 <- 0 

d2xprime . dade <- dxprime. de/a 

d2yprime.dade <- dyprime. de/a 

d2xprime . dadtheta <- dxprime . dtheta/a 

d2yprime . dadtheta <- dyprime . dtheta/a 

ddenom.de <- -2 * e * cossq. tt . theta 

ddenom. dtheta <- - e A 2 * sin(2 * (tt - theta)) 

terml <- ( - a A 2 * sinsq. 2 . tt . theta) /4 

xprime . denom. sq <- xprime * denom A 2 

d2xprime.de2 <- xprime . denom. sq - e * (dxprime.de * denom A 2 + 2 * 
xprime * denom * ddenom.de) 

d2xprime.de2 <- ( terml/xprime . denom. sq A 2 ) * d2xprime.de2 
d2xprime.de2 [abs (xprime) < le-006] <- 0 

d2yprime.de2 <- one . minus . e . sq * (-2 * xprime * d2xprime.de2 - 2 * 

dxprime . de A 2 ) + 8 * e * xprime * dxprime.de - 2 * (a A 2 - 
xprime A 2 ) 

terml <- - a A 2 * e 

d2xprime . dedtheta <- - xprime . denom. sq * sin(4 * (tt - theta)) - 

sinsq . 2 . tt . theta * (xprime * denom * ddenom. dtheta + 
dxprime . dtheta * denom A 2) 

d2xprime . dedtheta <- (( - a A 2 * e ) /xprime . denom. sq A 2 ) * 
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( 



d2xpnme . dedtheta 

d2xprime . dedtheta [abs (xprime) < le-006] <- 0 

d2yprime . dedtheta <- (-l/yprime A 2 ) * (yprime * ( (one .minus . e . sq * 



### 

### 

### 



xprime * d2xprime . dedtheta + dxprime . dtheta * dxprime.de)) - 
2 * e * xprime * dxprime . dtheta ) - one .minus . e . sq * xprime * 
dxprime . dtheta * dxprime.de) 
d2yprime . dedtheta [abs (yprime) < le-006] <- 0 
dnum. dtheta <- -2 * one . minus . e . sq * cos (2 * (tt - theta)) 
d2xprime . dtheta2 <- xprime . denom. sq * dnum. dtheta - num * (2 * 

xprime * denom * ddenom. dtheta + denom A 2 * dxprime . dtheta) 
d2xprime . dtheta2 <- (d2xprime . dtheta2 * a A 2)/(2 * 

xprime . denom. sq A 2 ) d2xprime. dtheta2 [abs (xprime) < le-006] <-0 
d2yprime . dtheta2 <- - ( one .minus . e . sq/ yprime ) * (yprime * (xprime 

* d2xprime . dtheta2 + dxprime . dtheta A 2 ) - dyprime . dtheta * ( 

xprime * dxprime . dtheta) ) 
d2yprime . dtheta2 [abs (yprime) < le-006] <- 0 

d2x.dade <- cos. theta * d2xprime . dade - sin. theta * d2yprime . dade 
d2y.dade <- sin. theta * d2xprime . dade + cos.thet 
d2x.de2 <- cos. theta * d2xprime.de2 - sin. theta 
d2y.de2 <- sin. theta * d2xprime.de2 + cos. theta 
grad.a2 <- 2 * sum(dx.da A 2 + dy.da A 2) 

grad.ae <- 2 * sum( - x.diff * d2x.dade + dx.da * dx.de - y.diff 
d2y.dade + dy.da * dy.de) 

d2x.dadtheta <- cos. theta * d2xprime . dadtheta - sin. theta * 
d2yprime . dadtheta - dy.da 

d2y. dadtheta <- sin. theta * d2xprime . dadtheta + cos. theta * 
d2yprime . dadtheta + dx.da 

d2xprime . dedtheta - sin. theta * 



* d2yp rime . dade 
d2 yprime . de2 
d2 yprime. de2 



d2x. dedtheta <- cos. theta 

d2yprime . dedtheta - dy.de 
d2y. dedtheta <- sin. theta * d2xprime . dedtheta 
d2yprime . dedtheta + dx.de 
d2x.dtheta2 <- cos. theta 
d2yprime . dtheta2 - 
d2y.dtheta2 <- sin. theta 
d2yprime . dtheta2 + 



+ cos. theta * 



2 

★ 

2 



d2xprime . dtheta2 - sin. theta * 

* dy. dtheta + x 

d2xprime . dtheta2 + cos. theta * 

* dx. dtheta + y 
grad.atheta <- 2 * sum( - x.diff * d2x. dadtheta + dx.da * 
dx. dtheta - y.diff * d2y. dadtheta + dy.da * dy. dtheta) 
grad.e2 <- 2 * sum ( - x.diff * d2x.de2 + dx.de A 2 - y.diff * 

d2y.de2 + dy.de A 2) 

grad.etheta <- 2 * sum( - x.diff * d2x. dedtheta + dx.de * 
dx. dtheta - y.diff * d2y. dedtheta + dy.de * dy. dtheta) 
grad.theta2 <- 2 * sum( - x.diff * d2x.dtheta2 + dx.dtheta A 2 
y.diff * d2y.dtheta2 + dy.dtheta A 2) # 
grad. mat <- matrix (c (grad. a2, grad.ae, grad.atheta, grad.ae, 
grad.e2, grad.etheta, grad.atheta, grad.etheta, 
grad. theta2 ) , 3, 3, T) 
cat (’’Hessian approx, is. . .\n") 
print (solve (grad. mat) ) 
if ( fit . center == F) 

hessian <- c(grad.a2, grad.ae, grad.e2, grad.atheta, 
grad. e theta, grad. theta2 ) 

else { 



Second derivatives: a and xO, a and yO 

d2xprime . dadxO <- dxprime . dxO/a 
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=#: 



d2 yp rime . dadxO <- dyprime . dxO/ a 
d2 xp rime . dadyO <- dxprime . dyO/a 
d2yprime . dadyO <- dyprime . dyO/a 

d2x. dadxO <- cos. theta * d2 xp rime . dadxO - sin. theta * 
d2yprime . dadxO 

d2y. dadxO <- sin. theta * d2xprime . dadxO + cos. theta * 
d2 yp rime . dadxO # 

# and the gradient 

grad. axO <- 2 * sum( - x.diff * d2x. dadxO + dx.da * dx.dxO 
y.diff * d2y. dadxO + dy.da * dy.dxO) 
d2x. dadyO <- cos. theta * d2 xp rime . dadyO - sin. theta * 
d2yprime . dadyO 

d2y. dadyO <- sin. theta * d2xp rime . dadyO + cos. theta * 
d2 yp rime . dadyO # 

# 

grad.ayO <- 2 * sum( - x.diff * d2x. dadyO dx.da * dx.dyO 
y.diff * d2y. dadyO + dy.da * dy.dyO) # 

# 

# 

ddenom.dxO <- e A 2 * sin . 2 . tt . theta * dt.dxO 
ddenom.dyO <- e A 2 * sin . 2 . tt . theta * dt.dyO 
A <- xprime * denom A 2 

dA.dxO <- denom * (2 * xprime * ddenom.dxO + denom * 

dxprime . dxO) 

dA.dyO <- denom * (2 * xprime * ddenom.dyO + denom * 

dxprime . dyO ) 

out. front <- - (a A 2 * e)/4 

z <- sinsq. 2 . tt . theta/A # 
z[abs(A) < tol] <- 0 

d2 xprime . dedxO <- out. front * ( (2 * sin(4 * (tt - theta)) * 

dt.dxO) /A - ((dA.dxO * z ) /A) ) 

d2xprime . dedyO <- out. front * ( (2 * sin(4 * (tt - theta)) * 

dt.dyO) /A - ((dA.dyO * z)/A)) 
d2xprime . dedxO [ abs (A) < tol] <- 0 
d2xprime . dedyO [abs (A) < tol] <- 0 # 

# Here's one from Mathematica. 

# 

d2yp rime. dedxO O ( - (3 * e) /2 * yprime * dt.dxO * 
sin. 2 . tt . theta) /denom A 2 
d2 yp rime . dedxO [abs (yprime) < tol] <- 0 
d2x. dedxO <- cos. theta * d2xp rime . dedxO - sin. theta * 
d2 yprime . dedxO 

d2y. dedxO <- sin. theta * d2 xp rime . dedxO + cos . theta * 
d2 yprime . dedxO 

grad.exO <- 2 * sum( - x.diff * d2x. dedxO + dx.de * dx.dxO 
y.diff * d2y. dedxO + dy.de * dy.dxO) 
d2 yp rime . dedyO <- ( - (3*e)/2* yprime * dt.dyO * 

sin . 2 . tt . theta ) /denom A 2 
d2 yp rime . dedyO [abs (yprime) < tol] <- 0 
d2x. dedyO <- cos. theta * d2 xp rime . dedyO - sin. theta * 
d2 yp rime . dedyO 

d2y. dedyO <- sin. theta * d2xprime . dedyO + cos. theta * 
d2 yprime . dedyO 

grad.eyO <- 2 * sum( - x.diff * d2x. dedyO + dx.de * dx.dyO 
y.diff * d2y. dedyO + dy.de * dy.dyO) # 

Here's another from Mathematica. 
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# 

# 



# 

# 



# 



# 



out. front <- (xprime * (1 - 2 * e A 2 + e A 2 * 

cos . 2 . tt . theta) ) / denom A 2 
out . front [abs (denom) < tol] <- 0 
d2xprime . dthetadxO <- out. front * dt . dxO 
d2xprime.dthetady0 <- out. front * dt.dyO # 

out. front <- ( (one .minus . e . sq) * yprime * (3 - 2 * denom))/ 

denom A 2 

out . front [abs (denom) < tol] <- 0 
d2yprime . dthetadxO <- out. front * dt.dxO 
d2yprime . dthetadyO <- out. front * dt.dyO # 



d2x. dthetadxO <- - dy.dxO + cos. theta * d2xprime . dthetadxO 

- sin. theta * d2yprime . dthetadxO 
d2y . dthetadxO <- (dx.dxO - 1) + sin. theta * 

d2xprime. dthetadxO + cos. theta * d2yprime . dthetadxO 
grad.thetaxO <- 2 * sum( - x.diff * d2x . dthetadxO + 
dx.dtheta * dx.dxO - y.diff * d2y . dthetadxO + 
dy.dtheta * dy.dxO) 

d2x . dthetadyO <- - (dy.dyO - 1) + cos. theta * 

d2xprime. dthetadyO - sin. theta * d2yp rime . dthetadyO 
d2y . dthetadyO <- dx.dyO + sin. theta * d2xprime . dthetadyO + 
cos. theta * d2yprime . dthetadyO 
grad.thetayO <- 2 * sum( - x.diff * d2x. dthetadyO + 
dx.dtheta * dx.dyO - y.diff * d2y. dthetadyO + 
dy.dtheta * dy.dyO) 

# 



d2t . dx02 <- -2 * (dt.dxO * dt.dyO) 

d2t . dxOdyO < 1/R.sq + 2 * (dt.dxO) A 2 

d2t . dy02 <- - d2t.dx02 

d2xprime . dx02 <- - dxprime . dtheta * d2t.dx02 - 

d2xprime. dthetadxO * dt.dxO 
d2yprime . dx02 <- - dyprime . dtheta * d2t.dx02 - 

d2yprime. dthetadxO * dt.dxO 
d2x.dx02 <- cos. theta * d2xprime . dx02 - sin. theta * 
d2 yprime. dx02 

d2y.dx02 <- sin. theta * d2xprime . dx02 + cos. theta * 
d2yprime . dx02 

grad.x02 <- 2 * sum( - x.diff * d2x.dx02 + dx.dxO A 2 - y.diff 
* d2y . dx02 + dy.dxO A 2) # 



d2xprime. dxOdyO <- - dxprime . dtheta * d2t. dxOdyO - 

d2xprime . dthetadyO * dt.dxO 
d2 yprime. dxOdyO <- - dyprime . dtheta * d2t. dxOdyO - 

d2yprime . dthetadyO * dt.dxO 

d2x. dxOdyO <- cos. theta * d2xprime . dxOdyO - sin. theta * 
d2 yprime . dxOdyO 

d2y. dxOdyO <- sin. theta * d2xprime . dxOdyO + cos. theta * 
d2yprime . dxOdyO 

grad.xOyO <- 2 * sum( - x.diff * d2x. dxOdyO + dx.dxO * 
dx.dyO - y.diff * d2y. dxOdyO + dy.dxO * dy.dyO) # 



d2xprime . dy02 <- - dxprime . dtheta * d2t.dy02 - 

d2xprime. dthetadyO * dt.dyO 
d2yprime.dy02 <- - dyprime . dtheta * d2t.dy02 - 

d2yprime. dthetadyO * dt.dyO 
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d2x.dy02 <- cos. theta * d2xprime . dy02 - sin. theta * 
d2yprime . dy02 

d2y.dy02 <- sin. theta * d2xprime . dy02 + cos. theta * 
d2yprime . dy02 

grad.y02 <- 2 * sum( - x.diff * d2x.dy02 + dx.dyO A 2 - y.diff 
* d2y . dy02 + dy.dyO A 2) 

hessian <- c(grad.a2, grad.ae, grad.e2, grad.atheta, 
grad.etheta, grad.theta2, grad.axO, grad.exO, 
grad. thetaxO, grad.x02, grad.ayO, grad.eyO, 
grad. the tayO, grad.xOyO, grad.y02) # 

### print (hessian) 

} 

thing <- list ( gradient = grad, hessian = hessian) 
return (thing) 



> ell . grad. II 

function (params, tt, my.x, my.y, is . there . hess, fit. center, class. I, 
which . type ) 

{ 

tol <- le-006 # 

a <- params [1] 
b <- params [2] 
e <- sqrt(l - (b/a) A 2) 
theta <- params [3] 
if ( fit .center == T) { 

center. x <- params [4] 
center. y <- params [5] 
delta <- params [6] 

tt <- ell. tt (my.x - center. x, my.y - center. y) 

} 

else { 

delta <- params [4] 
center. x <- center. y <- 0 

} 

if (sum (class . I ) == length (my . x) ) { 

fitted <- ell.pred(tt, a, b, theta, return . unrotated. too = 
T, fit. center = fit. center, center. x = center. x, center. y = 
center . y ) 

xprime <- fitted$x . prime 
yprime <- fitted$y. prime 
x <- fitted$x 
y <- fitted$y 

} 

else { 

fitted. I <- ell . pred ( tt [class . I] , a, b, theta, 

return . unrotated. too = T, fit. center = fit. center, 
center. x = center. x, center. y = center. y) 
if (which. type == 1) 

fitted. II <- ell.pred(tt [ Iclass.I] , a + delta, b, 
theta, return. unrotated. too = T, fit. center 
= fit. center, center. x = center. x, center. y 
= center. y) 

else fitted. II <- ell . pred ( tt [! class . I ] , a, b + delta, 
theta, 

return . unrotated. too = T, fit. center = 
fit. center, center. x = center. x, center. y = 
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center . y) 

xprime <- yprime <- x <- y <- numeric (length (my . x) ) 

xprime [class . I] <- fitted . I$x . prime 

xprime [ ! class . I] <- fitted. II$x. prime 

yprime [class . I] <- fitted. I$y. prime 

yprime [! class . I] <- fitted. II$y. prime 

x[class.I] <- fitted. I$x 

x[!class.I] <- fitted. II$x 

y[class.I] <- fitted. I$y 

y[!class.I] <- fitted. II$y 

} 

if (which. type == 1) { 

a <- rep(a, length (my . x) ) 
a[!class.I] <- a[!class.I] + delta 

} 

else { 

b <- rep(b, length (my . x) ) 
b[!class.I] <- b[!class.I] + delta 

} 

cos. theta <- cos (theta) 
sin. theta <- sin (theta) 

cos . 2 . tt . theta <- cos(2 * (tt - theta)) 
sin. 2 . tt . theta <- sin(2 * (tt - theta)) 
sinsq. tt . theta <- (sin(tt - theta) ) A 2 
cossq. tt. theta <- (cos(tt - theta) ) A 2 
sinsq. 2 . tt . theta <- (sin(2 * (tt - theta) )) A 2 
consq.2. tt.theta <- (cos (2 * (tt - theta) )) A 2 
one .minus . e . sq <- 1 - e A 2 

denom <- cossq. tt . theta * one .minus . e . sq + sinsq . tt . theta 
dxprime.da <- xprime/a 

dxprime.de <- ( ( - a A 2/(4 * xprime)) * (e * 
sinsq. 2 . tt . theta) ) / (denom A 2 ) 
dxp rime . de [abs (xprime) < tol] <- 0 
dyprime.da <- yprime/a 

dyprime.de <- - ( one .minus . e . sq * xprime * dxprime.de + e * (a A 2 

- xprime A 2) ) /yprime 
dyprime . de [abs (yprime) < tol] <- 0 # 



Okay. Here’s where we go from 



to 



We still need dx/y.de. 



de.db <- - b/(a A 2 

dx.da <- cos. theta 
dx.de <- cos. theta 

dx. db <- cos. theta 

de . db 

dy. da <- sin. theta 
dy.de <- sin. theta 
dy.db <- sin. theta 

de . db 

x. diff <- my.x - x 

y. diff <- my.y - y 
grad.a.indiv <- -2 
grad. a <- sum(grad. 
### grad.e <- -2 
grad.b.indiv <- -2 
grad.b <- sum(grad. 
if (which. type == 1) 

grad. delta <- 



* e) 

* dxprime.da - sin. theta * dyprime.da # 

* dxprime.de - sin. theta * dyprime.de 

* dxprime.de * de.db - sin. theta * dyprime.de * 

* dxprime.da + cos. theta * dyprime.da # 

* dxprime.de + cos. theta * dyprime.de 

* dxprime.de * de.db + cos. theta * dyprime.de * 



* (x.diff * dx.da + y.diff * dy.da) 

a . indiv) 

* sum (x.diff * dx.de + y.diff * dy.de) 

* (x.diff * dx.db + y.diff * dy.db) 

b . indiv) 

sum (grad. a . indiv [ ! class . I] ) 
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# 



else grad. delta <- sum ( grad. b . indiv[ ! class . I ] ) 
num <- one .minus . e . sq * sin(2 * (tt - theta}) 
dxprime.dtheta <- (a A 2/ (2 * xprime) ) * (num/denom A 2) 

dxprime. dtheta [abs (xprime) < tol] <- 0 

dyp rime . dthet a <- - (one .minus . e . sq * xprime * dxprime.dtheta)/ 

yprime 

dyp r ime. dtheta [abs (yp rime) < tol] <- 0 

dx. dtheta <- - (y - center. y) + cos. theta * dxprime.dtheta - 

sin. theta * dyp rime . dtheta 

dy. dtheta <- (x - center. x) + sin. theta * dxprime.dtheta + 
cos. theta * dyprime . dtheta 

grad. theta <- -2 * sum(x.diff * dx. dtheta + y.diff * dy. dtheta) 



# if (! exists ( "killme" , frame = 0)) 

browser ( ) 

if ( fit . center == F) 

grad <- c (grad. a, grad.b, grad. theta, 

else { 

dxprime. dt <- - dxprime.dtheta 

dyprime. dt <- - dyprime . dtheta 

R. sq <- (my.x - center. x) A 2 + (my.y - 
dt.dxO <- (my.y - center . y) /R. sq 
dt.dyO <- - (my.x - center . x) /R . sq 

dxprime. dxO <- dxprime. dt * dt.dxO 
dxprime. dyO <- dxprime. dt * dt.dyO 
dyp rime. dxO <- dyprime. dt * dt.dxO 
dyprime. dyO <- dyprime. dt * dt.dyO # 
dyprime .dxprime <- one .minus . e . sq * ( 

dyprime .dxprime [abs (yprime) < tol] <- 

dx. dxO <- (cos. theta * dxprime. dxO) - 

dyprime. dxO) + 1 

dy. dxO <- (sin. theta * dxprime. dxO) + 

dyprime .dxO) 

dx. dyO <- (cos. theta * dxprime. dyO) - 

dyprime. dyO) 

dy. dyO <- (sin. theta * dxprime. dyO) + 

dyprime. dyO) + 1 

grad.xO <- -2 * sum(x.diff * dx.dxO + 
grad.yO <- -2 * sum(x.diff * dx.dyO + 
grad <- c (grad. a, grad.b, grad. theta, 
grad. delta ) 

} 



grad. delta) 



center .y) A 2 



- xprime/ yprime) 

0 # 

(sin. theta * 

(cos. theta * 

(sin. theta * 

(cos. theta * 

y.diff * dy.dxO) 
y.diff * dy.dyO) 
grad.xO, grad.yO, 



if ( is . there . hess == F) { 
print (grad) 
return (grad) 

} 

d2xprime.da2 <- d2yprime.da2 <- 0 
d2xprime.dade <- dxprime. de/a 
d2 yprime. dade <- dyp rime.de/ a 
d2xprime. dadtheta <- dxprime . dtheta/a 
d2yprime . dadtheta <- dyprime . dtheta/a 
ddenom.de <- -2 * e * cossq. tt . theta 
ddenom. dtheta <- - e A 2 * sin(2 * (tt - theta)) 

terml <- ( - a A 2 * sinsq. 2 . tt . theta) /4 
xprime . denom. sq <- xprime * denom A 2 

d2xprime.de2 <- xprime . denom. sq - e * (dxprime.de * denom A 2 + 
xprime * denom * ddenom.de) 

d2xprime.de2 <- ( terml/xprime . denom. sq A 2 ) * d2xprime.de2 



2 * 
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### 

### 

### 

### 



### 

### 



d2xpr±me .de2 [abs (xprime) < le-006] <- 0 

d2yprime . de2 <- one .minus . e . sq * (-2 * xprime * d2xprime.de2 - 2 * 

dxprime . de A 2 ) + 8 * e * xprime * dxprime.de - 2 * (a A 2 - 
xprime A 2 ) 

terml <- - a A 2 * e 

d2xprime . dedtheta <- - xprime . denom. sq * sin(4 * (tt - theta)) - 

sinsq. 2 . tt . theta * (xprime * denom * ddenom. dtheta + 
dxprime . dtheta * denom A 2 ) 

d2xprime . dedtheta <- ( ( - a A 2 * e ) /xprime . denom. sq A 2 ) * 

d2 xprime . dedtheta 

d2xprime. dedtheta [abs (xprime) < le-006] <- 0 

d2yp rime . dedtheta <- (-l/yprime A 2 ) * (yprime * ( (one .minus . e . sq * 

xprime * d2xprime . dedtheta + dxprime . dtheta * dxprime.de)) - 
2 * e * xprime * dxprime . dtheta ) - one .minus . e . sq * xprime * 

dxprime . dtheta * dxprime.de) 
d2yprime . dedtheta [abs (yprime) < le-006] <- 0 
dnum. dtheta <- -2 * one .minus . e . sq * cos (2 * (tt - theta)) 
d2xprime . dtheta2 <- xprime . denom. sq * dnum. dtheta - num * (2 * 

xprime * denom * ddenom. dtheta + denom A 2 * dxprime . dtheta ) 
d2xprime .dtheta2 <- (d2xprime . dtheta2 * a A 2)/(2 * 
xprime . denom. sq A 2 ) 

d2xprime.dtheta2 [abs (xprime) < le-006] <- 0 

d2yprime . dtheta2 <- - (one .minus . e . sq/ yprime) * (yprime * (xprime 

* d2xprime . dtheta2 + dxprime . dtheta A 2 ) - dyprime . dtheta * ( 

xprime * dxprime . dtheta ) ) 

d2yprime.dtheta2 [abs (yprime) < le-006] <- 0 # 

d2e . db2 <- (b * de . db - e)/(a * e) A 2 
d2e . dadb <- (2 * b) / (a A 3 * e) # 

d2x.dade <- cos. theta * d2xprime . dade - sin. theta * d2yprime . dade 
d2x.dadb <- cos. theta * (dxprime.de * d2e.dadb + d2xprime . dade * 
de.db) - sin. theta * (dyprime.de * d2e.dadb + d2yprime . dade 

* de.db) # 

d2y.dade <- sin. theta * d2xprime . dade + cos. theta * d2yprime . dade 
d2y.dadb <- sin. theta * (dxprime.de * d2e.dadb + d2xprime . dade * 
de.db) + cos. theta * (dyprime.de * d2e.dadb + d2yprime . dade 

* de.db) # 

d2x.de2 <- cos. theta * d2xprime . de2 - sin. theta * d2yprime.de2 
d2y.de2 <- sin. theta * d2xprime.de2 + cos. theta * d2yprime.de2 
d2x.dadb <- cos. theta * (dxprime.de * d2e.dadb + d2xprime . dade * 
de.db) - sin. theta * (dyprime.de * d2e.dadb + d2yprime . dade 

* de.db) 

d2y.dadb <- sin. theta * (dxprime.de * d2e.dadb + d2xprime . dade * 
de.db) + cos. theta * (dyprime.de * d2e.dadb + d2yprime . dade 

* de.db) 

d2xprime . dedb <- d2xprime.de2 * de.db 
d2yprime . dedb <- d2yprime.de2 * de.db 

d2x.db2 <- cos. theta * (dxprime.de * d2e.db2 + d2xprime . dedb * 

de.db) - sin. theta * (dyprime.de * d2e.db2 + d2yprime . dedb * 
de.db) 

d2y.db2 <- sin. theta * (dxprime.de * d2e.db2 + d2xprime . dedb * 

de.db) + cos. theta * (dyprime.de * d2e.db2 + d2yprime . dedb * 
de. db) 

grad.a2 <- 2 * sum(dx.da A 2 + dy.da A 2) # 

grad.ae <- 2 * sum( - x.diff * d2x.dade + dx.da * dx.de - y.diff * 
d2y.dade + dy.da * dy.de) 

grad.ab <- 2 * sum( - x.diff * d2x.dadb + dx.da * dx.db - y.diff * 
d2y.dadb + dy.da * dy.db) 
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sin. theta * 



d2x.dadtheta <- cos. theta 
d2yprime . dadtheta 
d2y.dadtheta <- sin. theta 



d2xprime . dedtheta 
* de.db - dy.db 



de.db - sin. theta 



de.db + cos. theta 



d2xprime . dadtheta - 
dy . da 

d2xp rime. dadtheta + cos. theta * 
d2yprime. dadtheta + dx.da # 

### d2x. dedtheta <- cos. theta * d2xp rime . dedtheta - sin. theta * 

### d2yp rime . dedtheta - dy.de 

### d2y. dedtheta <- sin. theta * d2xprime . dedtheta + cos . theta * 

### d2yp rime . dedtheta + dx.de 

d2x.dbdtheta <- cos. theta 

* d2 yprime . dedtheta 

d2y.dbdtheta <- sin. theta * d2xp rime . dedtheta 

* d2 yprime . dedtheta * de.db + dx.db 
d2x.dtheta2 <- cos. theta * d2xprime . dtheta2 - sin. theta * 

d2yprime . dtheta2 - 2 * dy.dtheta + x 
d2y.dtheta2 <- sin. theta * d2xprime . dtheta2 + cos. theta * 
d2yprime.dtheta2 + 2 * dx.dtheta + y 
grad . atheta <- 2 * sum( - x.diff * d2x. dadtheta + dx.da * 
dx.dtheta - y.diff * d2y. dadtheta + dy.da * dy.dtheta) 

### grad.e2 <- 2 * sum( - x.diff * d2x.de2 + dx.de A 2 - y.diff * ### 
d2y . de2 + dy.de A 2) 

### grad.etheta <- 2 * sum( - x.diff * d2x. dedtheta + dx.de *dx.dtheta 
### - y.diff * d2y. dedtheta + dy.de * dy.dtheta) 

grad . b2 <- 2 * sum( - x.diff * d2x.db2 + dx.db A 2 - y.diff * 
d2y . db2 + dy.db A 2) 

grad.btheta <- 2 * sum( - x.diff * d2x.dbdtheta + dx.db * 

dx.dtheta - y.diff * d2y.dbdtheta + dy.db * dy.dtheta) 
grad.theta2 <- 2 * sum( - x.diff * d2x.dtheta2 + dx.dtheta A 2 - 
y.diff * d2y.dtheta2 + dy.dtheta A 2) 
if ( fit . center == F) 

hessian <- c(grad.a2, grad.ab, grad.b2, grad, atheta, 
grad.btheta, grad.theta2) 

else { 



# 



Second derivatives: a and xO, a and yO 

d2xprime . dadxO <- dxprime . dxO/a 
d2yprime . dadxO <- dyprime . dxO/a 
d2xprime . dadyO <- dxprime . dyO/a 
d2yprime.dady0 <- dyprime . dyO/a 

d2x. dadxO <- cos. theta * d2xp rime . dadxO - sin. theta * 
d2 yprime . dadxO 

d2y. dadxO <- sin. theta * d2xp rime . dadxO + cos . theta * 
d2yprime . dadxO # 

# and the gradient 

grad.axO <- 2 * sum ( - x.diff * d2x. dadxO + dx.da * dx.dxO - 
y.diff * d2y. dadxO + dy.da * dy.dxO) 
d2x. dadyO <- cos . theta * d2xp rime . dadyO - sin. theta * 
d2yprime . dadyO 

d2y. dadyO <- sin. theta * d2xp rime . dadyO + cos. theta * 
d2yprime . dadyO # 

# 

grad.ayO <- 2 * sum( - x.diff * d2x. dadyO + dx.da * dx.dyO - 
y.diff * d2y. dadyO + dy.da * dy.dyO) # 

# 

# 

ddenom.dxO <- e A 2 * sin. 2 . tt . theta * dt.dxO 
ddenom.dyO <- e A 2 * sin . 2 . tt . theta * dt.dyO 
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# Here's 

# 



### 

### 

### 

### 

### 

### 



### 

### 

### 

### 



# 

# Here's 

# 



# 



A <- xprime * denom A 2 

dA.dxO <- denom * (2 * xprime * ddenom.dxO + denom * 

dxprime.dxO) 

dA.dyO <- denom * (2 * xprime * ddenom.dyO + denom * 
dxprime.dyO) 

out. front <- - (a A 2 * e)/4 

2 O sinsq. 2 . tt . theta/A # 
z[abs(A) < tol] <- 0 

d2xprime . dedxO <- out. front * ( (2 * sin (4 * (tt - theta)) * 

dt . dxO ) /A - ((dA.dxO * z)/A)) 

d2xprime.dedy0 <- out. front * ( (2 * sin (4 * (tt - theta)) * 

dt.dyO)/A - ((dA.dyO * z)/A)) 
d2 xprime . dedxO [abs (A) < tol] <- 0 
d2xprime . dedyO [abs (A) < tol] <- 0 # 

one from Mathematica. 



d2 yp rime . dedxO <- ( - (3 * e) /2 
sin. 2 . tt. theta) /denom A 2 
d2yprime. dedxO [abs (yp rime) < tol] <- 0 
d2x. dedxO <- cos. theta * d2 xp rime . dedxO - 
d2yprime . dedxO 
d2y. dedxO <- sin. theta 
d2 yp rime . dedxO 

grad.exO <- 2 * sum( - x.diff 
y.diff * d2y. dedxO 
d2x.dbdx0 <- cos. theta * 

* d2 yp rime . dedxO * 
d2y.dbdx0 <- sin. theta * 

* d2yp rime. dedxO * 
grad.bxO <- 2 

y.diff * d2y . dbdxO 
d2 yp rime . dedyO <- ( - (3 

sin . 2 . tt . theta) /denom A 2 
d2yprime . dedyO [abs (yp rime) < tol] <- 0 
d2x. dedyO <- cos. theta * d2xprime . dedyO 
d2 ypr ime . dedyO 
d2y. dedyO <- sin. theta 
d2yprime . dedyO 
d2x.dbdy0 <- cos. theta 

* d2yprime . dedyO 
d2y.dbdy0 <- sin. theta 

* d2yp rime . dedyO 



yprime * dt.dxO * 



d2xprime . dedxO + 

d2x. dedxO 
+ dy.de * dy.dxO) 
d2 xp rime . dedxO * 
de . db 

d2xprime . dedxO * 
de . db 

sum ( - x.diff * d2x. dbdxO 
+ dy.db * dy.dxO) 
* e)/2 * yprime * 



d2 xprime . dedyO + 

d2xp rime . dedyO * 
de . db 

d2xprime . dedyO * 
de . db 



grad.byO <- 2 * sum( - x.diff * d2x.dbdy0 
y.diff * d2y.dbdy0 + dy.db * dy.dyO) 



sin. theta * 
cos. theta * 

+ dx.de * dx.dxO 
de.db - sin. theta 
de.db + cos. theta 
+ dx.db * dx.dxO 
dt.dyO * 

# 

sin. theta * 

cos. theta * 

de.db - sin. theta 

de.db + cos. theta 

+ dx.db * dx.dyO 

# 



another from Mathematica. 



out. front <- (xprime * (1 - 2 * e A 2 + e A 2 * 

cos . 2 . tt . theta) ) / denom A 2 
out . front [abs (denom) < tol] <- 0 
d2xprime.dthetadx0 <- out. front * dt.dxO 
d2xprime . dthetadyO <- out. front * dt.dyO # 

out. front <- ( (one .minus . e . sq) * yprime * (3 - 2 * denom))/ 

denom A 2 

out . front [abs (denom) < tol] <- 0 
d2yprime . dthetadxO <- out. front * dt.dxO 
d2yprime. dthetadyO <- out. front * dt.dyO # 
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d2x. dthetadxO <- - dy.dxO + cos. theta * d2xprime . dthetadxO 

- sin. theta * d2yprime . dthetadxO 
d2y. dthetadxO <- (dx.dxO - 1) + sin. theta * 

d2xprime. dthetadxO + cos. theta * d2 yp rime . dthetadxO 
grad . thetaxO <- 2 * sum( - x.diff * d2x . dthetadxO + 
dx.dtheta * dx.dxO - y.diff * d2y . dthetadxO + 
dy.dtheta * dy.dxO) 

d2x. dthetadyO <- - (dy.dyO - 1) + cos. theta * 

d2xprime. dthetadyO - sin. theta * d2yprime . dthetadyO 
d2y. dthetadyO <- dx.dyO + sin. theta * d2xprime . dthetadyO + 
cos. theta * d2yprime . dthetadyO 
grad.thetayO <- 2 * sum{ - x.diff * d2x . dthetadyO + 
dx.dtheta * dx.dyO - y.diff * d2y. dthetadyO + 
dy.dtheta * dy.dyO) 

# 

# 

# 

d2t . dx02 <- -2 * (dt.dxO * dt.dyO) 
d2t . dxOdyO <- -1/R.sq + 2 * (dt.dxO) A 2 
d2t . dy02 <- - d2t.dx02 

d2xprime . dx02 <- - dxprime . dtheta * d2t.dx02 - 

d2xprime. dthetadxO * dt.dxO 
d2yprime .dx02 <- - dyp rime . dtheta * d2t.dx02 - 

d2yprime. dthetadxO * dt.dxO 
d2x.dx02 <- cos. theta * d2xprime . dx02 - sin. theta * 
d2yprime. dx02 

d2y.dx02 <- sin. theta * d2xprime . dx02 + cos. theta * 
d2yprime. dx02 

grad.x02 <- 2 * sum( - x.diff * d2x.dx02 + dx.dxO A 2 - y.diff 
*d2y . dx02 + dy.dxO A 2) # 

# 

d2xprime. dxOdyO <- - dxprime . dtheta * d2t. dxOdyO - 

d2xprime. dthetadyO * dt.dxO 
d2yprime. dxOdyO <- - dyp rime . dtheta * d2t. dxOdyO - 

d2yprime . dthetadyO * dt.dxO 

d2x. dxOdyO <- cos. theta * d2xprime . dxOdyO - sin. theta * 
d2yprime . dxOdyO 

d2y. dxOdyO <- sin. theta * d2xp rime . dxOdyO + cos. theta * 
d2yprime . dxOdyO 

grad.xOyO <- 2 * sum( - x.diff * d2x. dxOdyO + dx.dxO * 
dx.dyO - y.diff * d2y. dxOdyO + dy.dxO * dy.dyO) # 

# 

d2xprime. dy02 <- - dxprime . dtheta * d2t.dy02 - 

d2xprime . dthetadyO * dt.dyO 
d2yprime. dy02 <- - dyprime . dtheta * d2t.dy02 - 

d2yprime . dthetadyO * dt.dyO 
d2x.dy02 <- cos. theta * d2xprime . dy02 - sin. theta * 
d2yprime.dy02 

d2y.dy02 <- sin. theta * d2xprime . dy02 + cos. theta * 
d2yprime . dy02 

grad.y02 <- 2 * sum( - x.diff * d2x.dy02 + dx.dyO A 2 - y.diff 
* d2y.dy02 + dy.dyO A 2) 

hessian <- c(grad.a2, grad.ab, grad.b2, grad.atheta, 
grad.btheta, grad.theta2, grad.axO, grad.bxO, 
grad. thetaxO, grad.x02, grad.ayO, grad.byO, 
grad.thetayO, grad.xOyO, grad.y02) # 

### hessian <- c(grad.a2, grad.ae, grad.e2, grad.atheta. 
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### grad.etheta, grad.theta2, grad.axO, grad.exO, 

### grad. thetaxO, grad.x02, grad.ayO, grad.eyO, 

### grad. thetayO, grad.xOyO, grad.y02) # 

### print (hessian) 

} 

thing <- list ( gradient = grad, hessian = hessian) 
return (thing) 

} 



> ell.pred 

function (tt, a, b, theta = 0, return . unrotated . too = F, fit. center = F, 
center. x = 0, center. y = 0) 

{ 

# 

# Get fitted x and y for ellipse with points at angles "tt", 

# with eccentricity = e and a = a. Rotate by theta afterwards, 

# if asked. Finally, if fit. center = T, move everything by 

# center. x in the x direction and by center. y in the y direction. 

# 

# A little algebra shows that 



# a A 2 sin A 2 (pi/2 - tt) (1 - e A 2) 

# x A 2 = if a > b. 

# sin A 2 ( tt ) + sin A 2 (pi/ 2 - tt) (1 - e A 2) 

# 

# (If a < b, that's y A 2, except you have to switch the tt ' s and the 

# (pi/2 - tt)'s.). The sin A 2 (pi/2-tt) term is "thang." So take 

# x (if a > b) to be the positive square root of that for the 



# moment. Then y A 2 = (a - ex) A 2 - (ae - x) A 2. So get that, too. 

# 

new.tt <- tt - theta 
if (a > b) { 

e <- sqrt ( 1 - (b/a) A 2) # 

### if(e > 0.99) return(1000 * length(x)) 

thang <- (sin(pi/2 - new.tt) A 2) * (1 - e A 2) # 

x <- sqrt((a A 2 * thang) /( sin (new. tt ) A 2 + thang)) 
yy <- (a - e * x) A 2 - (a * e - x) A 2 # 

# Make sure y A 2 is always positive (round-off errors can hurt here) ; 

# then get y. 

# 

yy [yy < o] <- o 

y <- sqrt (yy) # 

} 

else { 

e <- sqrt ( 1 - (a/b) A 2) # 

### if(e > 0.99) return(1000 * length(x)) 

thang <- (sin (new. tt ) A 2 ) * (1 - e A 2) # 

y <- b * sqrt (thang/ (sin (pi/2 - new.tt) A 2 + thang)) 
xx <- (b - e * y) A 2 - (b * e - y) A 2 # 
xx [xx < 0] <- 0 
x <- sqrt(xx) # 

) 

quad <- new.tt %% (2 * pi) 

quad. 2. 3 <- quad > pi/2 & quad < (3 * pi)/2 
x [ quad .2.3] <- - x[quad.2.3] 

quad. 3 . 4 <- quad > pi 
y[quad.3.4] <- - y[ quad. 3. 4] 

rotated. data <- matrix ( c ( cos ( theta ) , - sin (theta), sin (theta), 
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3*: 



cos ( theta)), 2, 2 , T) %*% rbind(x, y) 
if ( fit . center == F) { 
center. x <- 0 
center. y <- 0 

} 

if (return . unrotated. too == T) 

return (list (x = rotated. data [ 1 , ] + center. x, y = 

rotated. data [ 2 , ] + center. y, x. prime = x, y. prime 

= y) ) 

else return (list (x = rotated. data [ 1, ] + center. x, y = 

rotated. data [2, ] + center. y) ) 

} 



> ell.tt 

function (x, y) 

{ 



ell.tt: get angle values for data. It's acos (x/r) 
in the first quadrant, etc. 

tt <- numeric (length (x) ) 

ratio <- x/sqrt(x A 2 + y^2) 

ratio [ratio >1] <- 1 

ratio [ratio < -1] <- -1 

ind <- x >= 0 & y >= 0 

tt[ind] <- acos ( ratio [ind] ) 

ind <- x < 0 & y >= 0 

tt[ind] <- pi - acos ( - ratio [ind]) 

ind Ox<0&y<0 

ttfind] <- pi + acos ( - ratio [ind]) 

ind <- x >= 0 & y < 0 

tt[ind] <- 2 * pi - acos ( ratio [ind] ) # 

tt 
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APPENDIX B. EXAMPLES 



pclssl = pc is subject l’s initials 

the 1 is the spatial frequency 

ssl means sum of squares obtained from ellipse (vice ellipse.II or 
ellipse.ni) 

pclall = a data frame containing all of the data at 1 cpd for pc 

this data always consists of non-oblique data followed by oblique data. 
For this subject the first 80 (x,y) pairs are non-oblique and 81-160 are 
oblique data 

fit.center=T lets the ellipse center “float” vice being pinned to the origin 

> pc 1 ss 1 _ellipse(pc 1 all [, 1 ], pclall[,2], fit.center=T) 

a: 0.0615 ,b: 0.05194 ,th: -0.005236 ;x,y: 0.001551 0.00008819 ;obj: 0.06455 
a: 0.02742 ,b: 0.02742 ,th: -0.07295 ;x,y: 0.01894 -0.001528 ;obj: 0.07366 
removed middle output to save space 

a: 0.04698 ,b: 0.03466 ,th: -0.02604 ;x,y: 0.004495 0.0006379 ;obj: 0.02068 

a: 0.04698 ,b: 0.03466 ,th: -0.02605 ;x,y: 0.004495 0.0006379 ;obj: 0.02068 

a: 0.04698 ,b: 0.03466 ,th: -0.02605 ;x,y: 0.004495 0.0006379 ;obj: 0.02068 

> pclssl 
Message was 

RELATIVE FUNCTION CONVERGENCE 
a b theta center.x center.y 
0.04698378 0.03466318 -0.02604936 0.004495335 0.0006378881 

> pclssl[2] 

Sobjective: 

[1] 0.02068425 

pclss2.m = pc is subject l’s initials 

the 1 is the spatial frequency 

ss2 means sum of squares obtained from ellipse.II or ellipse.DI 
•III lets us know this is from ellipse.IQ 

> pclss2.III_ellipse.III (pclall[,l], pclall[,2], grad=F, 

+ is. there. hess=F, fit.center=T, class.I = (1:160) <81, plot=F) 
a: 0.07282 , delta.a: 0 ,b: 0.0615 delta.b: 0 ,th: -0.005236 ;x,y: 0.001551 
0.00008819 ;obj: 0.1364 

a: 0.07283 , delta.a: 0 ,b: 0.0615 delta.b: 0 ,th: -0.005236 ;x,y: 0.001551 
0.00008819 ;obj: 0.1364 

removed middle output to save space 
6 ;x,y: 0.004467 0.0006382 ;obj: 0.02067 

a: 0.04703 , delta.a: -0.0001285 ,b: 0.03424 delta.b: 0.0008527 ,th: -0.0266 
6 ;x,y: 0.004467 0.0006382 ;obj: 0.02067 
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a: 0.04703 , delta.a: -0.0001281 ,b: 0.03424 delta.b: 0.000853 ,th: -0.02666 
;x,y: 0.004467 0.0006382 ;obj: 0.02067 

a: 0.04703 , delta.a: -0.0001281 ,b: 0.03424 delta.b: 0.0008523 ,th: -0.0266 
6 ;x,y: 0.004467 0.0006382 ;obj: 0.02067 
Warning messages: 

1: singularity encountered in: nlminb.0(temp, p, liv, lv, objective, bounds, 
scale) 

removed identical warnings numbered 2 and 3 
4: singularity encountered in: nlminb.0(temp, p, liv, lv, objective, bounds, 
scale) 

> pclss2.III 
Message was 

RELATIVE FUNCTION CONVERGENCE 
a b theta center.x center.y delta.a 
0.04703421 0.03424435 -0.0266577 0.004466774 0.0006382473 -0.0001281251 
delta.b 

0.0008526696 

this is the objective function value which the program minimizes 

> pclss2.HI[2] 

Sobjective: 

[1] 0.02067104 

p-values for the ellipses being different? 

> 1 -pf(((pc 1 ss 1 [ [2] ] -pc 1 ss2.m[[2] ])/2)/(pc 1 ss 1 [[2]]/( 1 60-5)), 2, 155) 

[1] 0.9517269 

> l-pf(((pc3ssl[[2]]-pc3ss2.III[[2]])/2)/(pc3ssl[[2]]/(160-5)),2,155) 

[1] 0.07660905 

> 1 -pf(((pc7ss 1 [ [2] ]-pc7ss2. HI[[2] ])/2)/(pc7ss 1 [ [2 ] ]/( 1 60-5)), 2, 155) 

[1] 1.822829e-007 



are the b’s significant 

> l-pf((pclss2.n[[2]]-pclss2.III[[2]])/(pclss2.III[[2]]/l 53), 1,153) 
[1] 0.7693711 

> l-pf((pc3ss2.II[[2]]-pc3ss2.III[[2]])/(pc3ss2.III[[2]]/l 53), 1,153) 
[1] 0.837154 

> l-pf((pc7ss2.II[[2]]-pc7ss2.III[[2]])/(pc7ss2.ffl[[2]]/153), 1,153) 
[1] 2.110824e-008 



changing which. type to 2 forces the non-oblique and oblique ellipses to have the 
same major axis. 
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> pclss2.II.which.type.2_ellipse.II(pclall[,l],pclall[,2],fit.center=T,class.I = (1:160) < 
8 1 ,which.type=2,grad=F,is.there.hess=F,plot.it=F) 

>pc3 ss2.II.which.type.2_ellipse.II(pc3all[,l],pc3all[,2],fit.centen r T,class.I = (1:160) < 
81, which. type=2,grad=F, is. there. hess=F, plot. it=F) 

>pc7ss2.II.which.type.2_ellipse.II(pc7all[,l],pc7all[,2],fit.center=T,class.I = (1:160) < 
81,which.type=2,grad=F,is.there.hess=F,plot.it=F) 



are the a’s significant ? 
subject 1 at 7 cpd 

> l-pf((pc7ss2.H.which.type.2[[2]]-pc7ss2.ni[[2]])/(pc7ss2.ffl[[2]]/153), 1,153) 
[1] 0.5883988 

subject 1 at 3 cpd 

> l-pf((pc3ss2.II.which.type.2[[2]]-pc3 ss2.III[[2]])/(pc3ss2.III[[2]]/l 53), 1,153) 
[1] 0.02835718 

subject 1 at 1 cpd 

> l-pf((pclss2.n.which.type.2[[2]]-pclss2.ni[[2]])/(pclss2.ni[[2]]/153), 1,153) 
[1] 0.9753455 
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